
#define _CRT_SECURE_NO_WARNINGS
#include <stdlib.h> 
#include <stdio.h>
#include <string.h>
#include <math.h>
#include <time.h>
#include "nrutil.h"
#include "mjlib.h"

#define NAGE 19
#define NTYPE 2
#define NSIZE 2
#define NMAKE 7
#define NCAT 2
#define nj 532
#define njnew 28
#define njused 504
#define maxt 50
#define ng 9

//global variable declarations

//toggles for main program
int dbprint = 0; // set to 1 to create debugging files
int verbose = 0; // set to 1 to print outputs and warnings as the solver runs
int rj_expectations_static = 0; // 0 updates rj/pj as model runs (this is the preferred option for the full model), 1 sets a static baseline rj/pj as an expectation (simpler equilibrium problem)
int scrap_expectations_static = 0; // 0 uses actual (transition) scrap rates as model runs (this is the preferred option for the full model), 1 sets a static baseline for scrap expectations (simpler equilibrium problem)
int depreciation_mode; //toggle  1: immediate  2: asset values  3: scheduled at original purchase
int imperfectcomp; //toggle  0: perfect competition among new car producers,  1: imperfect competition
int dqdpusedeqm; //toggle  0: dqdp from hh problem, 1 from full used eqm model
double tauexpectations; //value in [0,1] is fraction of future declines in tau that are used when forming expectations over (otherwise myopic) future rental prices
double mpgexpectations; //value in [0,1] is fraction of future improvements in mpg that are used when forming expectations over (otherwise myopic) future rental prices

//total income, VMT and gas price (for the representative consumers in both regions)
double M_1, M_2; //current period income (unadjusted for profits)
double M_base_1, M_base_2; //initial period income
double M_growth; //exogenous income growth
double M_total_1, M_total_2;
double gas_price;

//data for cars
double mpg1[nj+1], mpg2[nj+1];
double basempg1[nj+1], basempg2[nj+1];
double commonmpg[nj+1]; //this is the fixed cost-related component of the fuel economy of a car

//indices for use accessing car vectors
int make[nj+1], age[nj+1], size[nj+1], type[nj+1]; //relate age, make, size and type to the car vector indexed j
int mc[NMAKE+1][NTYPE*NSIZE+1]; //links car vector to make and vehicle for the new car problem


//basic global vectors (prices, demands, etc.)
double rj1[nj+1], rj2[nj+1]; //rental prices
double exprj1[nj+1][maxt+NAGE+1], exprj2[nj+1][maxt+NAGE+1]; //expected future rental prices

double baserj1[nj+1], baserj2[nj+1]; //base case rental prices

double basepj1[nj+1], basepj2[nj+1]; //base purchase prices
double lastpj1[nj+1], lastpj2[nj+1]; //last step's purchase prices

double tapj1[nj+1],tapj2[nj+1]; //total annual operating cost

double d1[nj+1], d2[nj+1]; //demands
double d1_when_new[nj+1], d2_when_new[nj+1]; //original quantity when new
double based1[nj+1], based2[nj+1]; //base case demands
double lastd1[nj+1], lastd2[nj+1]; //last step's demands

double pscrap1[nj+1], pscrap2[nj+1]; //probability a car is scrapped beginning of current time period
double basepscrap1[nj+1], basepscrap2[nj+1]; //base case probability a car is scrapped beginning of current time period

double exppscrap1[nj+1], exppscrap2[nj+1]; //expected probability of scrappage at beginning of next time period
double lastexppscrap1[nj+1], lastexppscrap2[nj+1]; //expected probability of scrappage at beginning of current time period, as of end of previous time period
double exphj1[nj+1], exphj2[nj+1]; //expectations for hj
double lastexphj1[nj+1], lastexphj2[nj+1];
double exppj1[nj+1], exppj2[nj+1]; //expectations for pj
double lastexppj1[nj+1], lastexppj2[nj+1];

double hj1[nj+1], hj2[nj+1]; //expected repair cost paid for each not-scrapped vehicle (paid at beginning of period)
double basehj1[nj+1], basehj2[nj+1];

double payment1[nj+1][maxt+NAGE+1],payment2[nj+1][maxt+NAGE+1]; //stores the payment due on vehicle j in time t for certain versions of the depreciation calculation
double fracremain0[nj+1]; //baseline fraction of original sales remaining (assumes steady state)
double pj1[nj+1], pj2[nj+1]; //purchase prices
double valj1[nj+1], valj2[nj+1]; //value of (new) vehicles to used vehicle rental company

double maxsup1[nj+1], maxsup2[nj+1]; //maximum supply of a used vehicle (unscrapped ones from the previous period)
double usedsup1[nj+1], usedsup2[nj+1]; //for single used vehicle market, this keeps tracks of the actual supply of used vehicles of both types in the national market

double taxrev_tot_1, taxrev_tot_2; // circulation tax revenues

//vmt
double vmt[nj + 1]; // VMT by model-age schedule

//pollution parameters
double phi_co[nj + 1]; // CO emissions rates in grams per mile - current time step
double phi_hc[nj + 1]; // HC emissions rates in grams per mile - current time step
double phi_nox[nj + 1]; // NOx emissions rates in grams per mile - current time step
double phi_co_when_new[nj + 1]; // CO emissions when the vehicle was new - current time step
double phi_hc_when_new[nj + 1]; // HC emissions when the vehicle was new - current time step
double phi_nox_when_new[nj + 1]; // NOx emissions when the vehicle was new - current time step

//j x t matrix of pollution and parameters needed to construct it
int input_phi_maxt; //number of time steps to directly read phi (must read at least 1 time period to calibrate)
double jt_phi_co[nj + 1][maxt+NAGE+1]; // CO emissions rates in grams per mile, all j x t
double jt_phi_hc[nj + 1][maxt+NAGE+1]; // HC emissions rates in grams per mile, all j x t
double jt_phi_nox[nj + 1][maxt+NAGE+1]; // NOx emissions rates in grams per mile, all j x t
double jt_phi_co_when_new[nj + 1][maxt+NAGE+1]; // CO emissions when the vehicle was new, all j x t
double jt_phi_hc_when_new[nj + 1][maxt+NAGE+1]; // HC emissions when the vehicle was new, all j x t
double jt_phi_nox_when_new[nj + 1][maxt+NAGE+1]; // NOx emissions when the vehicle was new, all j x t

double phi_co_newcarstep[njnew + 1]; // CO emissions rates in grams per mile -- exogenous step for new cars
double phi_hc_newcarstep[njnew + 1]; // HC emissions rates in grams per mile -- exogenous step for new cars
double phi_nox_newcarstep[njnew + 1]; // NOx emissions rates in grams per mile -- exogenous step for new cars
double phi_co_agefactor[nj + 1]; // CO emissions age factor (multiplied by emissions when new to get current emissions)
double phi_hc_agefactor[nj + 1]; // HC emissions age factor (multiplied by emissions when new to get current emissions)
double phi_nox_agefactor[nj + 1]; // NOx emissions age factor (multiplied by emissions when new to get current emissions)


double damage_co, damage_hc, damage_nox; // pollution damages in dollars per metric ton
double damagepercar[nj + 1]; // pollution damages per vehicle (in a given time step) in dollars
double lcdamagepercar[nj + 1]; // pollution damages per vehicle (in a given time step), with production emissions attached to age 0 cars
double proddamage[nj + 1]; // production-associated damages per car (does not vary with time) in dollars per vehicle

// pollution aggregates
double co[nj + 1]; // CO emissions level in metric tons (aggregate, gets multiplied by Q)
double hc[nj + 1]; // HC emissions level in metric tons (aggregate, gets multiplied by Q)
double nox[nj + 1]; // NOx emissions level in metric tons (aggregate, gets multiplied by Q)
double totdamage[nj+1]; // total damage (aggregate, gets multiplied by Q)
double lctotdamage[nj+1]; // total damage above, with production emissions included

//circulation taxes
double tau[nj + 1];
double tau_read[nj+1][maxt+NAGE+1];

//pre-existing property tax
double proptaxrate; // pre-existing annual property tax rate: 0.01 = 1% of p
double proptax[nj + 1]; //tax per car in dollars


//vectors and arrays for the producer problem
//policy variables
int policytype; //0 no tau; 1 fractional Pigou; 2 fractional age-level Pigou; 3 new car only fractional lifetime Pigou; 4 flatten proptax revneutral; 5 additional proptax; 6 fractional Pigou rev neutral; 7 age cutoff rev nuetral; 8 age cutoff; 10 read from file
int policystart; //1 start in time step 1 = 2000, 9 start in time step 9=2016
int policyagestart; //age at which the binary age-based policy takes effect (level set using pollutionpolicy)
int lcpolicy; // 0 means policy is only on VMT damages, 1 is life cycle (attach production damage to age=0)
double pollutionpolicy; //tax set to a multiple of damages (e.g. 1 is for tax at 100% of damages)
double agepolicy; //tax set to a multiple of damages for a given age weighted by baseline demand (e.g. 1 is for tax at 100% of damages at age)
double newcarpolicy; //tax only on new cars, set to a fraction of discounted lifetime damages
double proptaxpolicy; //property tax change as a fraction of vehicle value, e.g. -0.015 is negative 1.5%
double proptaxconstant; // constant term in property tax assessment, adjusts to keep revenue neutral
double feebconstant; // constant term in feebate version, adjusts to keep revenue neutral
double basetailpipecost[njnew+1][maxt+NAGE+1]; //cost of base tailpipe standard in dollars
double tailpipecost[nj+1], jt_tailpipecost[njnew+1][maxt+NAGE+1]; //cost of tailpipe standard in dollars (total, so could just reflect baseline or have a policy cost added)
double tailpipez, tailpipev[njnew+1]; // tailpipe cost function parameters
double tailpipew[njnew+1][maxt+NAGE+1]; //baseline shifter in tailpipe function

double tailpipepolicy[njnew+1][maxt+NAGE+1]; // .01 means a 1% improvement, i.e. all emissions for that vintage (including as it ages) are reduced by 1%

double cafelimT[NMAKE+1], cafelimC[NMAKE+1], cafelimP[NMAKE+1]; //limits by make
double cafeC[NMAKE+1], cafeT[NMAKE+1], cafeP[NMAKE+1]; //actual average fuel efficiency by make
double cafeC_baset1[NMAKE+1], cafeT_baset1[NMAKE+1], cafeP_baset1[NMAKE+1]; //actual average fuel efficiency by make in base time 1
double cafeCstat[NMAKE+1], cafeTstat[NMAKE+1]; //set to 1 if CAFE binds for this car/truck & make, 0 if not
double cafePstat[NMAKE+1]; //1 if region-level standard is binding (not used in this model version)
int binding_violations; //variable that checks if no constraints are violated
double pavtarget[maxt+1]; //region-level target for each step (not used in this model version)
double cafeCtarget[maxt+1], cafeTtarget[maxt+1]; //uniform CAFE target for each step
int splitmpg = 0; //set to 1 to allow split mpg (in and out)
int splitused = 1; //set to 1 to solve for two separate used car markets (in and out)
//lagrange multipliers
double lambdaC[NMAKE+1], lambdaT[NMAKE+1], lambdaP[NMAKE+1]; //lagrange multipliers on the three constraints
double lambdafocs[NMAKE*3+1]; //stores lambda-related FOC's (i.e. value of the transformed constraints)
double lambdaCread[NMAKE+1], lambdaTread[NMAKE+1]; //calibrated values for lambda in baseline
//derivatives of demand and elasticities
double dqdp1[njnew+1][njnew+1], dqdp2[njnew+1][njnew+1]; //new car demand derivs
double **dqdrnew1, **dqdrnew2; //new car rental price derivs
double dqdp_elast1[njnew+1][njnew+1], dqdp_elast2[njnew+1][njnew+1]; //new car price elasticities
double dqdr_elast_all_1[nj+1][nj+1], dqdr_elast_all_2[nj+1][nj+1]; //all car rental price elasticities
double dqde1[nj+1][nj+1], dqde2[nj+1][nj+1]; //new car fuel economy derivs
//other
double cost1[nj+1], cost2[nj+1], basecost1[nj+1], basecost2[nj+1];  //calibrated cost variables
double costz[nj+1]; //fixed cost for common mpg technology
double ace[nj+1], bce[nj+1]; //cost of mpg technology.  cost = ace*x + bce*x^2  where x is improvement in mpg over baseline
double acz[nj+1], bcz[nj+1]; //cost of common mpg technology.
double argmincost1[nj+1], argmincost2[nj+1], argmincostz[nj+1]; //argmins of cost funcs
double mincost1[nj+1], mincost2[nj+1], mincostz[nj+1]; //mins of cost funcs
double dcde1[nj+1], dcde2[nj+1]; //first derivative of cost function wrt mpg, at current fuel economy
double dcdz1[nj+1], dcdz2[nj+1]; //first derivative of cost function wrt common mpg, at current fuel economy
double dczdz[nj+1]; //first derivative of the common mpg cost function wrt common mpg
double profit1, baseprofit1, profit2, baseprofit2; //total profits from new car production (difference is added to income since base assumed already included)
double profitrental1, baseprofitrental1, profitrental2, baseprofitrental2; //total profits from car rental market (difference is added to income since base assumed already included)
double profit_adj1, profit_adj2; //adjustment to income due to changes in (new and used) profit relative to baseline
double proptaxrev1, proptaxrev2, baseproptaxrev1, baseproptaxrev2; //property tax revenue
double exogmpgtech; //annual percent exogenous right shift in cost function
double theta_1; //fraction of common technologies (for calibration)

//derived variables
double milesremaining1[nj+1], milesremaining2[nj+1]; //expected, discounted, miles remaining (see scrap_expectations) for each car
double dvalde1[nj+1], dvalde2[nj+1]; //derivative of value of car j with respect to mpg (new cars)

//elasticity (rho) parameters (same across regions)
double rho_1[NAGE + 1][NTYPE + 1][NSIZE + 1];
double rho_2[NTYPE + 1][NSIZE + 1];
double rho_3[NTYPE + 1];
double rho_4; 
double rho_5;

//ownership cost and (composite) total annual car prices (ownership + gasoline)
//NB: the last subscript in the name refers to the region 
double tap_1_1[NAGE + 1][NTYPE + 1][NSIZE + 1][NMAKE + 1]; //matrix version of rj1 below
double tap_1_2[NAGE + 1][NTYPE + 1][NSIZE + 1][NMAKE + 1]; //matrix version of rj2 below
double tap_comp_1_1[NAGE + 1][NTYPE + 1][NSIZE + 1];
double tap_comp_1_2[NAGE + 1][NTYPE + 1][NSIZE + 1];
double tap_comp_2_1[NTYPE + 1][NSIZE + 1];
double tap_comp_2_2[NTYPE + 1][NSIZE + 1];
double tap_comp_3_1[NTYPE + 1];
double tap_comp_3_2[NTYPE + 1];
double tap_comp_4_1; //refers to p_v
double tap_comp_4_2;
double tap_x_1; //refers to p_x
double tap_x_2;
double tap_comp_5_1; //refers to p*: the "overall" price index
double tap_comp_5_2;

//vehicle demands
double dem_1_1[NAGE + 1][NTYPE + 1][NSIZE + 1][NMAKE + 1];
double dem_1_2[NAGE + 1][NTYPE + 1][NSIZE + 1][NMAKE + 1];
double dem_comp_2_1[NAGE + 1][NTYPE + 1][NSIZE + 1];
double dem_comp_2_2[NAGE + 1][NTYPE + 1][NSIZE + 1];
double dem_comp_3_1[NTYPE + 1][NSIZE + 1];
double dem_comp_3_2[NTYPE + 1][NSIZE + 1];
double dem_comp_4_1[NTYPE + 1];
double dem_comp_4_2[NTYPE + 1];
double dem_comp_5_1[NCAT + 1];
double dem_comp_5_2[NCAT + 1];

//vehicle demand ratios
double ratio_1_1[NAGE + 1][NTYPE + 1][NSIZE + 1][NMAKE + 1];
double ratio_1_2[NAGE + 1][NTYPE + 1][NSIZE + 1][NMAKE + 1];
double ratio_2_1[NAGE + 1][NTYPE + 1][NSIZE + 1];
double ratio_2_2[NAGE + 1][NTYPE + 1][NSIZE + 1];
double ratio_3_1[NTYPE + 1][NSIZE + 1];
double ratio_3_2[NTYPE + 1][NSIZE + 1];
double ratio_4_1[NTYPE + 1];
double ratio_4_2[NTYPE + 1];

//distribution (alpha) parameters
double alpha_1_1[NAGE + 1][NTYPE + 1][NSIZE + 1][NMAKE + 1];
double alpha_1_2[NAGE + 1][NTYPE + 1][NSIZE + 1][NMAKE + 1];
double alpha_2_1[NAGE + 1][NTYPE + 1][NSIZE + 1];
double alpha_2_2[NAGE + 1][NTYPE + 1][NSIZE + 1];
double alpha_3_1[NTYPE + 1][NSIZE + 1];
double alpha_3_2[NTYPE + 1][NSIZE + 1];
double alpha_4_1[NTYPE + 1];
double alpha_4_2[NTYPE + 1];
double alpha_5_1[NCAT + 1];
double alpha_5_2[NCAT + 1];

//scrap probability parameters (same across regions)
double b[nj+1]; //constant term in scrap function (calibrated to reproduce baseline data)
double eta[nj+1]; //elasticity parameter

//counters and basic model operation toggles
int t=0; //time step
long ranseed; //random seed
int runmoderead; //read-in value of runmode (used for time period 1+)
int runmode = 0; //toggle for run mode, always starts in base mode for time period 0
// 0=base:  saves "base case" utility value
// 1=policy:  starts at the base solution prices
int nt; //number of time steps
int retvalnew; //return value from new car algorithm
int retvalused; //return value from used car algorithm
int retvaldda; //return value from deaggregation
int boundsviolations = 0; //sum of all bounds violations
int fullprodprob; //toggle = 1 if we are including full used car solution while evaluating new car focs
double sumfoc;

//file I/O global vars
char runname[512]; //text name of run
char baserunname[512]; //text name of corresponding base run (when doing a policy case)
char datadir[512]; //directory to read data from
char cardatafile[512], lambdafile[512], cafetargetfile[512], elasfile[512], damagesfile[512], pollutionfilename[512], pollutionstepfilename[512], pollutionagefactorfilename[512], vmtfilename[512], policyfilename[512], tailpipecostfilename[512], tailpipeparamsfilename[512], tailpipepolicyfilename[512]; //file names for inputs
char outdir[512]; //directory to write to
char sbuf[512]; //string buffer (used as temp storage for strings)
FILE *mainoutput, *debug, *debug2;

//variables for the equivalent variation calculation
double util1[maxt+1][3], util2[maxt+1][3]; //indirect utility levels
double ev1[maxt+1], ev2[maxt+1]; //equivalent variations

// variables for the expansion of vehicle numbers across income groups
double dda_base[nj+1][ng+1];
double dda[nj+1][ng+1];
double dda_jtarg[nj+1], dda_gtarg[ng+1]; //target vectors

//constants, misc.
double drate;  //discount rate

//symbols to be used as extern by fdjac or broydn2
int usedcar_in_jacobian = 0; //switch turned on when inside a used car jacobian
int prodprob_in_jacobian = 0; //turned on when inside a producer problem jacobian
int prodprob_exit_jacobian = 0; // 1 means producer problem has just exited jacobian
int fdjac_counter = 0; //counts number ofjacobian calls, used to print out jacobian values
int reuse_jacobian = 0; //external (used in solution algorithms) -- when equal to 1 the used car solver attempts to solve using previous derivatives
double rsave[1500][1500];
double jaceps = 0.001; //sets the value of epsilon used in the jacobian

//convergence measures
double sumedglobal, sumfocglobal;

//function prototypes (local)
void advance(int t);
void readinput(void);
void initialization(void);
void definecar(void);
void definecalib(void);
void newcarcalib(int t);
void assign_bindingness(int t);
void check_binding_violations(void);
void ncarbasecafe(int njd, double *guesscc, double *focs);
void newcarjacobians(void);
void newcardems_usedeqm(int njnd, double *newcarpj, double *newcaraggd);
void newcarjacobians_usedeqm(void);
void newcardems1(int njnd, double *newcarrj1, double *newcaraggd1);
void newcardems2(int njnd, double *newcarrj2, double *newcaraggd2);
void allcarjacobians(void);
void allcardems1(int njnd, double *newcarrj1, double *newcaraggd1);
void allcardems2(int njnd, double *newcarrj2, double *newcaraggd2);
void hhproblem(double *hhrj1, double *hhrj2, double *hhd1, double *hhd2, double hhM_1, double hhM_2);
void usedcarexdems(int nmkt, double *depju, double *exdemu);
void usedcareqm(void);
void eqm(void);
void writeout(int t);
void prodprobfocs(int numfocs, double *pparam, double *focs);
void indirectutil(double v_1, double x_1, double v_2, double x_2);
void dda_objective(int nummvec, double *mvec, double *diffs);

//function prototypes (numerical recipes)
void broydn (double x[], int n, int *check, void (*vecfunc)(int, double [], double []));
void broydn2(double x[], int n, int *check, void (*vecfunc)(int, double [], double []));
void fdjac(int n, double x[], double fvec[], double **df, void (*vecfunc)(int, double [], double []));

//************************************************************************************
//  START MAIN PROGRAM   
//************************************************************************************

int main (int argc, const char * argv[]) {
	FILE *baserjfile;
    char rline[512]; //temporary string buffer for line at a time reading (lines under 512 chars)
    int j, jx, ax, xt;
    double tempw, tempwavg;
	
	debug = fopen("../../../../../temp/debug.csv", "w");
	debug2 = fopen("../../../../../temp/debug2.csv", "w");
	
	initialization();
	readinput();
	definecar();
	definecalib(); //includes calibration program for utility function and scrap function parameters
	
	//set maximum supply for time period 1
	for (j=njnew+1;j<=nj;j++) {
		maxsup1[j] = 1.0 / (1.0 - basepscrap1[j]) * based1[j];
		maxsup2[j] = 1.0 / (1.0 - basepscrap2[j]) * based2[j];
	}

	//calibrate the new car problem (or read in the parameters if a policy case)
	newcarcalib(0); //doesn't impose additional cafe limits, just assigns cafe_baset1 values (the baseline fleet average)
	allcarjacobians(); // compute demand jacobians at base
    
    runmode = runmoderead; // switch to policy mode if read in
    
    //open main output file:
    sprintf(sbuf,"%s%s.csv",outdir,runname);
    mainoutput = fopen(sbuf, "w");
    
	
	//************ TIME LOOP ************
	for (t=1;t<=nt;t++) {  //main time loop
		sprintf(sbuf,"\nSTARTING TIME %d",t);  timestamp(sbuf);
		//read in the base rental prices and scrap probabilities if in the policy case
		if (runmode != 0) {
			sprintf(sbuf,"./basecalib/baserj%d.txt",t);
			baserjfile = fopen(sbuf, "r");
			for (j=1;j<=nj;j++) {	
				fscanf(baserjfile, "%lf %lf %lf %lf %lf %lf %lf %lf", &baserj1[j], &baserj2[j], &basepj1[j], &basepj2[j], &basepscrap1[j], &basepscrap2[j], &based1[j], &based2[j]);
			}
            fscanf(baserjfile, "%lf %lf", &baseprofit1, &baseprofit2);
            fscanf(baserjfile, "%lf %lf", &baseprofitrental1, &baseprofitrental2);
            fscanf(baserjfile, "%lf %lf", &baseproptaxrev1, &baseproptaxrev2);
            fscanf(baserjfile, "%s", baserunname); trim(baserunname);
            fclose(baserjfile);
		}
		
		//calculate damages per car and circulation taxes
		for (j = 1; j <= nj; j++) {
			damagepercar[j] = vmt[j] * (damage_co * phi_co[j] + damage_hc * phi_hc[j] + damage_nox * phi_nox[j]) / 1000000;
            // this is in dollars per time step (note vmt is per 2 year time step)
            lcdamagepercar[j] = damagepercar[j] + proddamage[j];
            tau[j] = 0.0;
		}
        
        if (policytype == 1 & t >= policystart) {
            for (j = 1; j <= nj; j++) {
                if (lcpolicy==1) tau[j] = pollutionpolicy * lcdamagepercar[j];
                else tau[j] = pollutionpolicy * damagepercar[j];
            }
        }
                
        if (policytype == 2 & t >= policystart) {
            //compute weighted (base q) average damage
            for (ax=0; ax<=18; ax++) {
                tempwavg = 0.0;
                tempw = 0.0;
                for (j = ax*njnew+1; j <= (ax+1)*njnew; j++) {
                    if (lcpolicy==1) tempwavg +=lcdamagepercar[j] * based1[j];
                    else tempwavg +=damagepercar[j] * based1[j];
                    tempw += based1[j];
                }
                tempwavg = tempwavg / tempw;
                for (j = ax*njnew+1; j <= (ax+1)*njnew; j++) {
                    tau[j] = agepolicy*tempwavg;
                }
            }
        }
        
        if (policytype == 3 & t >= policystart) {
            for (j=1; j<=njnew; j++){
                if (lcpolicy==1) tau[j] += proddamage[j];
                for (ax=0; ax<=18; ax++) {
                    jx = j+ax*njnew;
                    tau[j] += newcarpolicy * vmt[jx] * fracremain0[jx] * (1 / pow(1.0+drate,ax)) *
                    (damage_co * jt_phi_co[jx][t+ax] + damage_hc * jt_phi_hc[jx][t+ax] + damage_nox * jt_phi_nox[jx][t+ax]) / 1000000;
                }
            }
        }
        
        if (policytype == 4 & t >= policystart) {
            tempw = 0.0;
            proptaxconstant = 0.0;  // starting guess for proptaxconstant based on baseline composition (gets updated as model runs to ensure revenue neutrality)
            for (j=1; j<=nj; j++) {
                proptaxconstant += (based1[j]+based2[j]) * 2* proptaxpolicy * basepj1[j] * 1000;
                tempw+=(based1[j]+based2[j]);
            }
            proptaxconstant = -proptaxconstant / tempw;
            for (j=1; j<=nj; j++) {
                tau[j] = proptaxconstant + 2* proptaxpolicy * basepj1[j] * 1000;
            }
        }
        
        if (policytype == 5 & t >= policystart) {
            for (j=1; j<=nj; j++) {
                tau[j] = 2 * proptaxpolicy * basepj1[j] * 1000;
            }
        }
        
        if (policytype == 6 & t >= policystart) {
            tempw = 0.0;
            feebconstant = 0.0;  // starting guess for feebconstant based on baseline composition (gets updated as model runs to ensure revenue neutrality)
            for (j=1; j<=nj; j++) {
                if (lcpolicy==1) feebconstant += (based1[j]+based2[j]) * pollutionpolicy * lcdamagepercar[j];
                else feebconstant += (based1[j]+based2[j]) * pollutionpolicy * damagepercar[j];
                tempw+=(based1[j]+based2[j]);
            }
            feebconstant = -feebconstant / tempw;
            
            for (j = 1; j <= nj; j++) {
                if (lcpolicy==1) tau[j] = feebconstant + pollutionpolicy * lcdamagepercar[j];
                else tau[j] = feebconstant + pollutionpolicy * damagepercar[j];
            }
        }
        
        
        if (policytype == 7 & t >= policystart) {
            tempw = 0.0;
            feebconstant = 0.0;  // starting guess for feebconstant based on baseline composition (gets updated as model runs to ensure revenue neutrality)
            for (j=1; j<=nj; j++) {
                if (age[j] >= policyagestart) feebconstant += (based1[j]+based2[j]) * pollutionpolicy;
                tempw += (based1[j]+based2[j]);
            }
            feebconstant = -feebconstant / tempw;
                     
            for (j = 1; j <= nj; j++) {
                if (age[j] >= policyagestart) tau[j] = feebconstant + pollutionpolicy;
                else tau[j] = 0;
            }
        }
        
        if (policytype == 8 & t >= policystart) {

            for (j = 1; j <= nj; j++) {
                if (age[j] >= policyagestart) tau[j] = pollutionpolicy;
                else tau[j] = 0;
            }
        }
        
        
        
        if (policytype == 10 & t >= policystart) {
            //tau is read in directly (note it is in $, not $000s)
            for (j = 1; j <= nj; j++) tau[j] = tau_read[j][t];
        }
        
		assign_bindingness(t);
			
		//solve the model for a set of equilibrium prices in this period
		
		binding_violations=1;
		while (binding_violations!=0) {
			eqm();
			check_binding_violations();
		}
		
		//in the base case run, store r, p, pscrap, d, profits
        if (runmode == 0) {
			for (j=1;j<=nj;j++) {
				baserj1[j] = rj1[j]; //base case so we have solved for baserj
				baserj2[j] = rj2[j];
				basepj1[j] = pj1[j];
				basepj2[j] = pj2[j];
                basepscrap1[j] = pscrap1[j];
                basepscrap2[j] = pscrap2[j];
                basehj1[j] = hj1[j];
                basehj2[j] = hj2[j];
				based1[j] = d1[j];
				based2[j] = d2[j];
			}
			sprintf(sbuf,"./basecalib/baserj%d.txt",t);
			baserjfile = fopen(sbuf, "w");
			for (j=1;j<=nj;j++) {fprintf(baserjfile,"%20.12lf %20.12lf %20.12lf %20.12lf %20.12lf %20.12lf %20.12lf %20.12lf\n", baserj1[j], baserj2[j], basepj1[j], basepj2[j], basepscrap1[j], basepscrap2[j], based1[j], based2[j]);}
            fprintf(baserjfile, "%20.12lf %20.12lf\n", baseprofit1, baseprofit2);
            fprintf(baserjfile, "%20.12lf %20.12lf\n", baseprofitrental1, baseprofitrental2);
            fprintf(baserjfile, "%20.12lf %20.12lf\n", baseproptaxrev1, baseproptaxrev2);
            fprintf(baserjfile, "%s\n", runname);
			fclose(baserjfile);
		}
		
		//in the policy run for a national used car market, keep track of the used cars of each type
		if (runmode == 1 && splitused == 0){
			for (j=njnew+1;j<=nj;j++) {
				usedsup1[j]=maxsup1[j]*(1.0 - pscrap1[j]);
				usedsup2[j]=maxsup2[j]*(1.0 - pscrap2[j]);
			}
		}
		if(verbose) timestamp("Time period solved.");
		
		//writes the results to a file		
		writeout(t);
		
		if(verbose) printf("\nEV1 %lf", ev1[t]);
        if(verbose) printf("\nEV2 %lf", ev2[t]);
        if(verbose) printf("\n util1[t=%d][runmode=%d] = %lf", t, runmode, util1[t][runmode]);
        if(verbose) printf("\n util2[t=%d][runmode=%d] = %lf", t, runmode, util2[t][runmode]);
        if (dbprint) fprintf(debug,"end time period %d\n",t);
        
		//step all the cars and quantities forward to the next time period
		if (t<nt) advance(t);
		
	} //end of main time loop
	
    
    fclose(mainoutput);
	fclose(debug);
	fclose(debug2);
	timestamp("Finished.");
	return 0;
}

void newcarcalib(int t) {
    //calibrate the producer problem (e.g. cost functions for improving mpg)
    //note mpg1 is assumed = mpg2 in the base case and for calibration
	double *costg, *focvec;
	double  tmpd1, tmpd2, tmpd3;
	int retval, j, m;
	char rline[512];
	FILE *calf;
	
    if(verbose) printf("\nEntering newcarcalib()");
	
	sprintf(sbuf,"./basecalib/newcarcal%d.csv",t);
    calf = fopen(sbuf, "w");  //open file for writing
		  
    //first we need to get dq/dp at the base
    if (dqdpusedeqm==0) newcarjacobians(); //places the jacobian dq/dp into the matrix dqdp
    else newcarjacobians_usedeqm();

    //set time step 1 CAFE standard based on benchmark weighted avg fuel economy
    for (m=1;m<=NMAKE;m++) {
        
        cafeC_baset1[m] = ( d1[mc[m][1]] + d1[mc[m][2]] + d2[mc[m][1]] + d2[mc[m][2]] ) /
        ( (1.0/mpg1[mc[m][1]]) * d1[mc[m][1]] +
         (1.0/mpg1[mc[m][2]]) * d1[mc[m][2]] +
         (1.0/mpg2[mc[m][1]]) * d2[mc[m][1]] +
         (1.0/mpg2[mc[m][2]]) * d2[mc[m][2]]   );
        
        cafeT_baset1[m]=( d1[mc[m][3]] + d1[mc[m][4]] + d2[mc[m][3]] + d2[mc[m][4]] ) /
        ( (1.0/mpg1[mc[m][3]]) * d1[mc[m][3]] +
         (1.0/mpg1[mc[m][4]]) * d1[mc[m][4]] +
         (1.0/mpg2[mc[m][3]]) * d2[mc[m][3]] +
         (1.0/mpg2[mc[m][4]]) * d2[mc[m][4]]   );
        
        cafeP_baset1[m] = ( d1[mc[m][1]] + d1[mc[m][2]] + d1[mc[m][3]] + d1[mc[m][4]] ) /
        ( (1.0/mpg1[mc[m][1]]) * d1[mc[m][1]] +
         (1.0/mpg1[mc[m][2]]) * d1[mc[m][2]] +
         (1.0/mpg1[mc[m][3]]) * d1[mc[m][3]] +
         (1.0/mpg1[mc[m][4]]) * d1[mc[m][4]]   );
    }
    
    //assign the baseline (calibrated) values of lambda
    for (m=1;m<=NMAKE;m++) {
        lambdaC[m] = lambdaCread[m];
        lambdaT[m] = lambdaTread[m];
        lambdaP[m] = 0.0;
    }

    costg = dvector(1,njnew*2);
    focvec = dvector(1,njnew*2);

    //initialize costz to 0: no common added technologies in time step 0
    for (j=1;j<=njnew;j++) costz[j] = 0.0;

    //make starting guesses for calibration of costs (costg[j]) and deriv cost wrt fe (costg[njnew+j])
    for (j=1;j<=njnew;j++) costg[j] = 10.0; //calibrating one cost vector (regions always symmetric at base)
    for (j=1;j<=njnew;j++) costg[njnew+j] = 0.1;
  
    jaceps = 0.001; //
    broydn(costg,(2*njnew),&retval,ncarbasecafe);
    
    ncarbasecafe(njnew, costg, focvec); //run once at maximum to check
    if(verbose) printf("\nNew car calib done - return value %d", retval);
    if (dbprint) fprintf(debug,"\nNew car calib done - return value %d \n", retval);

    
    //compute (base) producer profits for time 0 (assumed to already be included in income)
    baseprofit1 = 0.0;
    baseprofit2 = 0.0;
    for (j=1;j<=njnew;j++) {
        // no subtraction for costz because it is 0 in time 0 by definition
        baseprofit1 += (basepj1[j] - cost1[j]) * d1[j];
        baseprofit2 += (basepj2[j] - cost2[j]) * d2[j];
    }
    profit1 = baseprofit1;
    profit2 = baseprofit2;

    //if there is tech with regional spillovers (i.e. common tech z) then compute values for a_z and b_z, and adjust bce (since a fraction theta of the technology embedded in current cars is assumed common)
    if (theta_1 > 0 && imperfectcomp == 1) {
        for (j=1;j<=njnew;j++) {
            acz[j] = ace[j] * (d1[j] + d2[j]);
            bcz[j] = (bce[j] / (2.0 * theta_1)) * (d1[j] + d2[j]);
            bce[j] = bce[j] / (2.0 * (1.0 - theta_1));
        }
    }
    else {
        for (j=1;j<=njnew;j++) {
            acz[j] = 0.0;
            bcz[j] = 0.0;
        }
    }
    
    //print results to calibration file
    fprintf(calf, "make, lidx, mpg1, cost1, costz, pj1, pj2, rj1, rj2, d1, d2, dqdp1, dqdp2, elast1, elast2, ace, acz, bce, bcz");
    for (j=1;j<=njnew;j++) {
        fprintf(calf, "\n%d, %d, %20.12lf, %20.12lf, %20.12lf, %20.12lf, %20.12lf, %20.12lf, %20.12lf, %20.12lf, %20.12lf, %20.12lf, %20.12lf, %20.12lf, %20.12lf, %20.12lf, %20.12lf, %20.12lf, %20.12lf", make[j], (type[j]-1)*2 + size[j],
                mpg1[j],
                cost1[j],
                costz[j],
                basepj1[j],
                basepj2[j],
                baserj1[j],
                baserj2[j],
                d1[j],
                d2[j],
                dqdp1[j][j],
                dqdp2[j][j],
                dqdp_elast1[j][j],
                dqdp_elast2[j][j],
                ace[j],
                acz[j],
                bce[j],
                bcz[j]);
    }
    //and print the lambdas and cafe limits:
    fprintf(calf, "\nmake, lambdatruck, lambdacar, lambdapav, cafeT, cafeC, cafeP");
    for (m=1;m<=NMAKE;m++) {
        fprintf(calf,"\n%d, %20.12lf, %20.12lf, %20.12lf, %20.12lf, %20.12lf, %20.12lf",m,lambdaT[m],lambdaC[m],lambdaP[m],cafeT_baset1[m],cafeC_baset1[m],cafeP_baset1[m]);
        cafeC[m] = cafeC_baset1[m];
        cafeT[m] = cafeT_baset1[m];
        cafeP[m] = cafeP_baset1[m];
    }
    //and print base profit
    fprintf(calf, "\nbase profit reg1");
    fprintf(calf,"\n%20.12lf",baseprofit1);
    fprintf(calf, "\nbase profit reg2");
    fprintf(calf,"\n%20.12lf",baseprofit2);
    
    free_dvector(costg,1,njnew);
    free_dvector(focvec,1,njnew);

	
	//set basempg (used in producer problem foc's to find costs)
	for (j=1;j<=njnew;j++) {
		basempg1[j] = mpg1[j];
		basempg2[j] = mpg2[j];
		basecost1[j] = cost1[j];
		basecost2[j] = cost2[j];
	}

	//find argmins of cost functions
	for (j=1;j<=njnew;j++) {
		argmincost1[j] = -1.0*ace[j] / (2.0 * bce[j]);
		mincost1[j] = basecost1[j] + ace[j]*(argmincost1[j]) + bce[j]*pow(argmincost1[j],2);
		argmincost2[j] = -1.0*ace[j] / (2.0 * bce[j]);
		mincost2[j] = basecost2[j] + ace[j]*(argmincost2[j]) + bce[j]*pow(argmincost2[j],2);
        if (theta_1>0 && imperfectcomp == 1) argmincostz[j] = -1.0*acz[j] / (2.0 * bcz[j]);
		if (theta_1>0 && imperfectcomp == 1) mincostz[j] = acz[j]*argmincostz[j] + bcz[j]*pow(argmincostz[j],2);
        fprintf(calf, "\nargmins for the cost functions");
		fprintf(calf, "\nargmin cost  %d  %lf",j,argmincost1[j]);
		fprintf(calf, "\nargmin costZ %d  %lf",j,argmincostz[j]);
	}

	fclose(calf);
}

void ncarbasecafe(int njd, double *guesscc, double *focs) {
	//takes values for cost and a_ce, returns value of FOCs for region 1 (for calibration)
	
	int j, k, l, lj;
	double sumfoc;
	double mC[njnew+1],mT[njnew+1],mP[njnew+1],dmC[njnew+1],dmT[njnew+1],dmP[njnew+1];
	double extrafocs[njnew+1];
	
	//put guesses in appropriate global spot
	for (j=1;j<=njnew;j++)  cost1[j] = guesscc[j];
	for (j=1;j<=njnew;j++)  cost2[j] = cost1[j];
	for (j=1;j<=njnew;j++)  ace[j] = guesscc[njnew+j];
	
	//compute dqde (discounted sum of expected value of an addl mile per gallon for all rental periods)

	for (j=1;j<=njnew;j++) {	
		for (k=1;k<=njnew;k++) {
			dqde1[j][k] = -dvalde1[k] * dqdp1[j][k];
			dqde2[j][k] = -dvalde2[k] * dqdp2[j][k];
		}
	}
	
	//calculate transformed fuel economies, uses region 1 since in calibration
	for (j=1;j<=njnew;j++) {
		
		//start with the "m" variables (and dm/de)
		if (type[j] == 1) {
			mC[j] = 1.0 - cafeC_baset1[make[j]] / mpg1[j];
			dmC[j] = cafeC_baset1[make[j]] * pow(mpg1[j],-2);
		}
		else			 { mC[j] = 0.0; dmC[j] = 0.0; }
		
		if (type[j] == 2) {
			mT[j] = 1.0 - cafeT_baset1[make[j]] / mpg1[j];
			dmT[j] = cafeT_baset1[make[j]] * pow(mpg1[j],-2);
		}
		else			 { mT[j] = 0.0; dmT[j] = 0.0; }
		
		mP[j] =  1.0 - cafeP_baset1[make[j]] / mpg1[j];
		dmP[j] = cafeP_baset1[make[j]] * pow(mpg1[j],-2);
	}
	
	//now compute the FOCs
    if (imperfectcomp == 1) {
        for (j=1;j<=njnew;j++) {
            
            //FOCs wrt price for region 1
            focs[j] = d1[j];
            extrafocs[j] = d2[j];
            for (l=1;l<=4;l++) {
                lj = mc[make[j]][l]; //lj is the index for vehicle 'l' made by firm 'make'
                focs[j] +=      ( pj1[lj] - cost1[lj] + lambdaC[make[lj]] * mC[lj] + lambdaT[make[lj]] * mT[lj] + lambdaP[make[lj]] * mP[lj] ) * dqdp1[lj][j];
                extrafocs[j] += ( pj2[lj] - cost2[lj] + lambdaC[make[lj]] * mC[lj] + lambdaT[make[lj]] * mT[lj] ) * dqdp2[lj][j];
            }
            
            //FOCs wrt fuel economy
            focs[njnew+j] = ( -ace[j] + lambdaC[make[j]] * dmC[j] + lambdaT[make[j]] * dmT[j] ) * (d1[j] + d2[j]) +
            lambdaP[make[j]] * dmP[j] * d1[j];
            for (l=1;l<=4;l++) {
                lj = mc[make[j]][l]; //lj is the index for vehicle 'l' made by firm 'make'
                focs[njnew+j] += ( pj1[lj] - cost1[lj] + lambdaC[make[lj]] * mC[lj] + lambdaT[make[lj]] * mT[lj] + lambdaP[make[lj]] * mP[lj] ) * dqde1[lj][j] +
                ( pj2[lj] - cost2[lj] + lambdaC[make[lj]] * mC[lj] + lambdaT[make[lj]] * mT[lj]  ) * dqde2[lj][j];
            }

        }
    }
    
    if (imperfectcomp == 0) {
        for (j=1;j<=njnew;j++) {
  
            //FOCs wrt price for region 1
            focs[j] = pj1[j] - cost1[j] + lambdaC[make[j]] * mC[j] + lambdaT[make[j]] * mT[j] + lambdaP[make[j]] * mP[j];
            extrafocs[j] =  pj2[j] - cost2[j] + lambdaC[make[j]] * mC[j] + lambdaT[make[j]] * mT[j] ;
            
            //FOCs wrt fuel economy
            focs[njnew+j] = ( -ace[j] + lambdaC[make[j]] * dmC[j] + lambdaT[make[j]] * dmT[j] + dvalde1[j]) * (d1[j] + d2[j]) + lambdaP[make[j]] * dmP[j] * d1[j];
        }
    }
	
	//scaling for solver
	for (j=1;j<=2*njnew;j++) {
		focs[j] = focs[j]*100000.0;
	}

	//debugging
	sumfoc = 0.0;
	for (j=1;j<=2*njnew;j++) sumfoc += fabs(focs[j]);
	sumfocglobal = sumfoc;
	if (usedcar_in_jacobian == 0) {
        if(verbose) printf("\nsumfoc calib %lf", sumfoc);
        if (dbprint) {
            for (j=1;j<=2*njnew;j++)  fprintf(debug, "%20.12lf, ", guesscc[j]);
            fprintf(debug,"\n");
            for (j=1;j<=2*njnew;j++)  fprintf(debug, "%20.12lf, ", focs[j]);
            for (j=1;j<=njnew;j++)  fprintf(debug, "%20.12lf, ", extrafocs[j]);
            fprintf(debug,"\n");
            for (j=1;j<=nj;j++)  fprintf(debug, "%20.12lf, ", d1[j]);
            fprintf(debug,"\n\n");
        }
	}
	
}


void newcarjacobians_usedeqm(void) {
	//writes the jacobian dqdp for new cars at current prices, solves used car eqm at each step to account for used vehicle price changes
	int j, k;
	double **dqdpall,*ptemp,*dtemp;
	static int ncjcount = 1;
	
    if(verbose) printf("\nEntering newcarjacob (with full eqm in used market)");
	
	dqdpall = dmatrix(1,2*njnew,1,2*njnew);
	ptemp=dvector(1,2*njnew);
	dtemp=dvector(1,2*njnew);
	
	jaceps = 0.005;

	for (j=1;j<=njnew;j++){
		ptemp[j]       = pj1[j];
		ptemp[njnew+j] = pj2[j];
	}

	fdjac(njnew*2, ptemp, dtemp, dqdpall, newcardems_usedeqm);

	for (j=1;j<=njnew;j++) {
		for (k=1;k<=njnew;k++) {
			
			dqdp1[j][k] = dqdpall[j][k];
			dqdp_elast1[j][k] = dqdpall[j][k] * ptemp[k] / dtemp[j]; 

			dqdp2[j][k] = dqdpall[njnew+j][njnew+k];
			dqdp_elast2[j][k] = dqdpall[njnew+j][njnew+k] * ptemp[njnew+k] / dtemp[njnew+j]; 

		}
	}

	ncjcount++;
		
	free_dmatrix(dqdpall,1,2*njnew,1,2*njnew);
	free_dvector(ptemp,1,2*njnew);
	free_dvector(dtemp,1,2*njnew);

	
}


void newcardems_usedeqm(int njnd, double *newcarpj, double *newcaraggd) {
	//takes a vector of new car purchase prices, solves for new car demands
	//used to provide dqdp to the producer problem
	int j;
    
    for (j=1;j<=njnew;j++) {
        pj1[j] = newcarpj[j];
        pj2[j] = newcarpj[njnew+j];
    }
    
	usedcareqm();

	for (j=1;j<=njnew;j++) {
		newcaraggd[j] = d1[j];
		newcaraggd[njnew+j] = d2[j];
	}
	
}



void newcarjacobians(void) {
	//compute the jacobian dqdp for new cars at current prices
	int j, k;
	double *tempdem1, *tempdem2;
	static int ncjcount = 1;
	
    if(verbose) printf("\nEntering newcarjacob()");
	
	tempdem1 = dvector(1,njnew);
	tempdem2 = dvector(1,njnew);
	
	jaceps = 0.001;

    //evaluate dqdr
	newcardems1(njnew,rj1,tempdem1); //get base new car demands for jacobian (at prices in rj1)
	fdjac(njnew,rj1,tempdem1,dqdrnew1,newcardems1);
    
	newcardems2(njnew,rj2,tempdem2); //get base new car demands for jacobian (at prices in rj2)
	fdjac(njnew,rj2,tempdem2,dqdrnew2,newcardems2);	

    //translate dqdr to dqdp
	for (j=1;j<=njnew;j++) {
		for (k=1;k<=njnew;k++) {
			if (rj_expectations_static==1) dqdp1[j][k] = dqdrnew1[j][k] * baserj1[k] / basepj1[k];
			else dqdp1[j][k] = dqdrnew1[j][k] * rj1[k] / pj1[k];
			dqdp_elast1[j][k] = dqdp1[j][k] * pj1[k] / tempdem1[j];
		}
	}
	for (j=1;j<=njnew;j++) {
		for (k=1;k<=njnew;k++) {
			if (rj_expectations_static==1) dqdp2[j][k] = dqdrnew2[j][k] * baserj2[k] / basepj2[k];
			else dqdp2[j][k] = dqdrnew2[j][k] * rj2[k] / pj2[k];
			dqdp_elast2[j][k] = dqdp2[j][k] * pj2[k] / tempdem2[j];   
		}
	}
		
	ncjcount++;
	free_dvector(tempdem1,1,njnew);	
	free_dvector(tempdem2,1,njnew);
	
}

void newcardems1(int njnd, double *newcarrj1, double *newcaraggd1) {
	//takes a vector of new car rental prices, solves for new car demands
	//used to provide dqdp to the producer problem
	int j;
	double rjpass1[nj+1], dempass1[nj+1];	
	
	for (j=1;j<=njnew;j++) rjpass1[j] = newcarrj1[j]; //guesses for new car rentals
	for (j=njnew+1;j<=nj;j++) rjpass1[j] = rj1[j]; //current (fixed) vector for used rentals
	
	hhproblem(rjpass1, rj2, dempass1, d2, M_1, M_2);   //get base demands

	for (j=1;j<=njnew;j++) newcaraggd1[j] = dempass1[j]; //just return new car demands, ignore the rest
	
}
void newcardems2(int njnd, double *newcarrj2, double *newcaraggd2) {
	//takes a vector of new car rental prices, solves for new car demands
	//used to provide dqdp to the producer problem
	int j;
	double rjpass2[nj+1], dempass2[nj+1];	
	
	for (j=1;j<=njnew;j++) rjpass2[j] = newcarrj2[j]; //guesses for new car rentals
	for (j=njnew+1;j<=nj;j++) rjpass2[j] = rj2[j]; //current vector for used rentals
	
	hhproblem(rj1, rjpass2, d1, dempass2, M_1, M_2);   //get base demands
	
	for (j=1;j<=njnew;j++) newcaraggd2[j] = dempass2[j]; //just return new car demands, ignore the rest
	
}

void allcarjacobians(void) {
	//compute the jacobians dqdr1 and dqdr2 for all (new and used) cars at current prices
	int j, k;
	double **dqdr1, **dqdr2;
	double *tempdem1, *tempdem2; //temp vector to store demands while in fdjac
	static int ncjcount = 1;
		
	if(verbose) printf("\nEntering allcarjacob()");
	
	dqdr1 = dmatrix(1,nj,1,nj);
	tempdem1 = dvector(1,nj);
	dqdr2 = dmatrix(1,nj,1,nj);
	tempdem2 = dvector(1,nj);
	
	jaceps = 0.001;

	allcardems1(nj,rj1,tempdem1); //get base car demands for jacobian (at prices rj1)
	fdjac(nj,rj1,tempdem1,dqdr1,allcardems1);

	allcardems2(nj,rj2,tempdem2); //get base car demands for jacobian (at prices rj2)
	fdjac(nj,rj2,tempdem2,dqdr2,allcardems2);

	for (j=1;j<=nj;j++) {
		for (k=1;k<=nj;k++) {
			dqdr_elast_all_1[j][k] = dqdr1[j][k] * rj1[k] / tempdem1[j];
		}
	}
	for (j=1;j<=nj;j++) {
		for (k=1;k<=nj;k++) {
			dqdr_elast_all_2[j][k] = dqdr2[j][k] * rj2[k] / tempdem2[j];   
		}
	}
	
	ncjcount++;
		
	free_dmatrix(dqdr1,1,nj,1,nj);
	free_dvector(tempdem1,1,nj);	
	free_dmatrix(dqdr2,1,nj,1,nj);
	free_dvector(tempdem2,1,nj);

}

void allcardems1(int njnd, double *allcarrj1, double *allcaraggd1) {
	//takes a vector of all car rental prices, returns all car demands
	//used to investigate all dqdr elasticities
	int j;
	double rjpass1[nj+1], dempass1[nj+1];	
	
	for (j=1;j<=nj;j++) rjpass1[j] = allcarrj1[j]; //guesses for all car rentals
	
	hhproblem(rjpass1, rj2, dempass1, d2, M_1, M_2);   //get demands
	
	for (j=1;j<=nj;j++) allcaraggd1[j] = dempass1[j];
	
}

void allcardems2(int njnd, double *allcarrj2, double *allcaraggd2) {
	//takes a vector of all car rental prices, returns all car demands
	//used to investigate all dqdp elasticities
	int j;
	double rjpass2[nj+1], dempass2[nj+1];	
	
	for (j=1;j<=nj;j++) rjpass2[j] = allcarrj2[j]; //guesses for all car rentals
		
	hhproblem(rj1, rjpass2, d1, dempass2, M_1, M_2);   //get demands
	
	for (j=1;j<=nj;j++) allcaraggd2[j] = dempass2[j];
	
}

void eqm(void) {
	//solves the FOCs for the producer problem (with the used vehicle market also clearing at each step)
	double *guess_n, *guess_nsave, *foc_n, ctC, ctT;
	int nmkt_n, j, m, x, maxitb, itb;
	
	nmkt_n = 4*njnew + 3*NMAKE; //solving for two sets of prices, one set of mpgs, one set of common mpgs, and 3 lambdas per make
	if (splitmpg == 1) nmkt_n += njnew; //add in the second set of mpgs if regional regulation
	guess_n = dvector(1,nmkt_n);
	guess_nsave = dvector(1,nmkt_n);
	foc_n = dvector(1,nmkt_n);
	
    //have solver compute a new jacobian (at the start of the time period)
    reuse_jacobian = 0;
    
    //trigger a restart of the derivative calculation
    fullprodprob=0;
	
	//make starting guesses for solution vector to new car problem
	for (j=1;j<=njnew;j++) {
		guess_n[j] = pj1[j]; //region 1 new car prices 
		guess_n[njnew+j] = pj2[j]; //region 2 new car prices
		guess_n[2*njnew+j] = mpg1[j];  //mpg1 (note that mpg2 just mirrors mpg1 if splitmpg=0)
		guess_n[3*njnew+j] = commonmpg[j];

	}
	for (m=1;m<=NMAKE;m++) {
		guess_n[4*njnew+m] = lambdaC[m];
		guess_n[4*njnew+m+NMAKE] = lambdaT[m];
		guess_n[4*njnew+m+NMAKE*2] = lambdaP[m];
	}	
	if (splitmpg == 1) {
		for (j=1;j<=njnew;j++) {
			guess_n[4*njnew+3*NMAKE+j] = mpg2[j];
		}
	}

	//evaluate the foc's once at the starting guesses, not solving them yet
    prodprobfocs(nmkt_n,guess_n,foc_n);
    if(verbose) printf("\nFirst evaluation of producer problem FOCs, intermediate boundsviolations %d",boundsviolations);

      
    jaceps = 0.005; // can adjust this to assist solver if model is not converging
    
    boundsviolations=0;
    broydn2(guess_n,nmkt_n,&retvalnew,prodprobfocs);
    
	if(verbose) printf("\nSolved producer problem FOC's retval %d intermediate boundsviolations %d",retvalnew,boundsviolations);

	//evaluate once more at solution (may be slightly updated with new profits and/or jacobian measure)
    boundsviolations=0;
	prodprobfocs(nmkt_n,guess_n,foc_n);

	printf("\n\nfull eqm retval %d final boundsviolations %d", retvalnew, boundsviolations);
	
	free_dvector(guess_n,1,nmkt_n);
	free_dvector(guess_nsave,1,nmkt_n);
	free_dvector(foc_n,1,nmkt_n);

	for (m=1;m<=NMAKE;m++) {
		cafeC[m] = ( d1[mc[m][1]] + d1[mc[m][2]] + d2[mc[m][1]] + d2[mc[m][2]] ) /
		( (1.0/mpg1[mc[m][1]]) * d1[mc[m][1]] +
		 (1.0/mpg1[mc[m][2]]) * d1[mc[m][2]] +
		 (1.0/mpg2[mc[m][1]]) * d2[mc[m][1]] +
		 (1.0/mpg2[mc[m][2]]) * d2[mc[m][2]]   );
		
		cafeT[m]=( d1[mc[m][3]] + d1[mc[m][4]] + d2[mc[m][3]] + d2[mc[m][4]] ) /
		( (1.0/mpg1[mc[m][3]]) * d1[mc[m][3]] +
		 (1.0/mpg1[mc[m][4]]) * d1[mc[m][4]] +
		 (1.0/mpg2[mc[m][3]]) * d2[mc[m][3]] +
		 (1.0/mpg2[mc[m][4]]) * d2[mc[m][4]]   );
		
		cafeP[m] = ( d1[mc[m][1]] + d1[mc[m][2]] + d1[mc[m][3]] + d1[mc[m][4]] ) /
		( (1.0/mpg1[mc[m][1]]) * d1[mc[m][1]] +
		 (1.0/mpg1[mc[m][2]]) * d1[mc[m][2]] +
		 (1.0/mpg1[mc[m][3]]) * d1[mc[m][3]] +
		 (1.0/mpg1[mc[m][4]]) * d1[mc[m][4]]   );
	}
}

void prodprobfocs(int numfocs, double *pparam, double *focs) {
	//takes guesses for new car prices, mpgs, and lambdas, returns value of all focs
	//note there are 4*njnew + 3*NMAKE foc's (2 sets of prices, one set of mpgs, one set of common mpgs, and 3 lambdas for each make)
	//if splitmpg=1, another njnew focs are added to the end for mpg2's
	int j, k, l, lj, m, n;
	double tmpdqdp1[njnew+1][njnew+1], tmpdqdp2[njnew+1][njnew+1], tmpd1[njnew+1], tmpd2[njnew+1];
	double mC1[njnew+1],mC2[njnew+1],mT1[njnew+1],mT2[njnew+1],mP1[njnew+1];
	double dmC1[njnew+1],dmC2[njnew+1],dmT1[njnew+1],dmT2[njnew+1],dmP1[njnew+1];	
	static int ppfcount = 0;
	FILE *tmp;

	// If evaluating producer problem jacobian (i.e. not yet searching for the solution), then can use demand approx.
	// but while actually solving, need the full solution to the used vehicle problem:
    if (prodprob_in_jacobian == 1) {
        fullprodprob = 0;
        prodprob_exit_jacobian = 0;
    }

    if (prodprob_in_jacobian == 0 & fullprodprob == 1) {
        prodprob_exit_jacobian = 0;
    }
    
    if (prodprob_in_jacobian == 0 & fullprodprob == 0) {
        fullprodprob=1;
        prodprob_exit_jacobian = 1;
    }
    
	//trap potential solver guesses of low/negative mpgs and prices
	
	for (j=1;j<=njnew;j++) {
		if (pparam[j] < 3.0) {  // PRICE 1
			printf("\nwarning: low or negative price guess for car1 %d", j);
			pparam[j] = 1.0 + 2.0 / exp(3.0-pparam[j]);
            boundsviolations++;
		}
		if (pparam[njnew+j] < 3.0) {
			printf("\nwarning: low or negative price guess for car2 %d", j);
			pparam[njnew+j] = 1.0 + 2.0 / exp(3.0-pparam[njnew+j]);
            boundsviolations++;
		}
		if (pparam[j+2*njnew] < 8.0) {  // MPG 1
			printf("\nMPG warning: low or negative mpg guess for car %d", j);
			pparam[j+2*njnew] = 6.0 + 2.0 / exp(8.0-pparam[j+2*njnew]);
            boundsviolations++;
		}
		if (splitmpg==1 && pparam[j+4*njnew + 3*NMAKE] < 8.0) {
			printf("\nMPG warning: low or negative mpg guess for car %d", j);
			pparam[j+4*njnew + 3*NMAKE] = 6.0 + 2.0 / exp(8.0-pparam[j+4*njnew + 3*NMAKE]);
            boundsviolations++;
		}
	}

    //copy guesses into pj (globals, needed in used car eqm)
    for (j=1;j<=njnew;j++) {
        pj1[j] = pparam[j];
        pj2[j] = pparam[njnew+j];
    }
    
	//place the guesses for new car fuel economy in global attribute array:
	for (j=1;j<=njnew;j++) {
		mpg1[j] = pparam[j + 2*njnew];
        if (theta_1 > 0 && imperfectcomp == 1) {
            commonmpg[j] = pparam[j + 3*njnew];
        }
        else {
            commonmpg[j] = 0.0;
        }
		if (splitmpg==1) mpg2[j] = pparam[j + 4*njnew + 3*NMAKE];
		else mpg2[j] = mpg1[j];
	}
    
    //finally take the pparams that represent lambdas and put them in the main lambda arrays
    for(m=1;m<=NMAKE;m++) {
        lambdaC[m] = cafeCstat[m] * pparam[4*njnew+m];
        lambdaT[m] = cafeTstat[m] * pparam[4*njnew+m+NMAKE];
        lambdaP[m] = cafePstat[m]  * pparam[4*njnew+m+NMAKE*2];
    }
    
	//calculate the derivative of the cost function at this level of fuel economy
	for (j=1;j<=njnew;j++) {
		dcde1[j] = ace[j] + bce[j]*2.0*(mpg1[j] - basempg1[j] - commonmpg[j]);
		dcde2[j] = ace[j] + bce[j]*2.0*(mpg2[j] - basempg2[j] - commonmpg[j]);

        
        if ( (mpg1[j]-basempg1[j]-commonmpg[j]) < argmincost1[j] ) {
            if(verbose) printf("\nCOST FUNCTION WARNING: searching below argmin1 for car %d", j);
            boundsviolations++;
        }
        if ( (mpg2[j]-basempg2[j]-commonmpg[j]) < argmincost2[j] ) {
            if(verbose) printf("\nCOST FUNCTION WARNING: searching below argmin2 for car %d", j);
            boundsviolations++;
        }
         
 
        if (theta_1 > 0 && imperfectcomp == 1) {
            dcdz1[j] = -1.0 * dcde1[j];
            dcdz2[j] = -1.0 * dcde2[j];
            dczdz[j] = acz[j] + bcz[j]*2.0*commonmpg[j];
            
            if ( commonmpg[j] < argmincostz[j] ) {
                if(verbose) printf("\nCOST FUNCTION WARNING: searching below argminz for car %d", j);
                boundsviolations++;
            }
        }
	}
	
    //calculate total cost for each vehicle
    for (j=1;j<=njnew;j++) {
        cost1[j] = tailpipecost[j]/1000 + basecost1[j] + ace[j]*(mpg1[j]-basempg1[j]-commonmpg[j]) + bce[j]*pow(mpg1[j]-basempg1[j]-commonmpg[j],2);
        cost2[j] = tailpipecost[j]/1000 + basecost2[j] + ace[j]*(mpg2[j]-basempg2[j]-commonmpg[j]) + bce[j]*pow(mpg2[j]-basempg2[j]-commonmpg[j],2);
    }

    //calculate total fixed cost
    if (theta_1 > 0 && imperfectcomp == 1) {
        for (j=1;j<=njnew;j++) {
            costz[j] = acz[j]*commonmpg[j] + bcz[j]*pow(commonmpg[j],2);
        }
    }
    else {
        costz[j] = 0.0;
    }
    
    //calculate fuel economies
    for (j=1;j<=njnew;j++) {
        if (type[j] == 1) {
            mC1[j] = 1.0 - cafelimC[make[j]] / mpg1[j];
            dmC1[j] = cafelimC[make[j]] * pow(mpg1[j],-2);
            mC2[j] = 1.0 - cafelimC[make[j]] / mpg2[j];
            dmC2[j] = cafelimC[make[j]] * pow(mpg2[j],-2);
        }
        else  {
            mC1[j] = 0.0;
            dmC1[j] = 0.0;
            mC2[j] = 0.0;
            dmC2[j] = 0.0;
        }
        
        if (type[j] == 2) {
            mT1[j] = 1.0 - cafelimT[make[j]] / mpg1[j];
            dmT1[j] = cafelimT[make[j]] * pow(mpg1[j],-2);
            mT2[j] = 1.0 - cafelimT[make[j]] / mpg2[j];
            dmT2[j] = cafelimT[make[j]] * pow(mpg2[j],-2);
        }
        else  {
            mT1[j] = 0.0;
            dmT1[j] = 0.0;
            mT2[j] = 0.0;
            dmT2[j] = 0.0;
        }
        
        mP1[j] =  1.0 - cafelimP[make[j]] / mpg1[j];
        dmP1[j] = cafelimP[make[j]] * pow(mpg1[j],-2);
    }
    
    //if perfect competition, price new cars at cost (so the solver is choosing only mpg):
    if (imperfectcomp == 0) {
        for (j=1;j<=njnew;j++) {
            pj1[j] = cost1[j] - lambdaC[make[j]] * mC1[j] - lambdaT[make[j]] * mT1[j] - lambdaP[make[j]] * mP1[j];
            pj2[j] = cost2[j] - lambdaC[make[j]] * mC2[j] - lambdaT[make[j]] * mT2[j];
        }
    }
        
    //now solve the used car problem at this set of new car prices and mpgs so that demands will update (we need FOCs to hold at eqm demands):
    usedcareqm();
    

    if (rj_expectations_static==0 & prodprob_exit_jacobian==1 & imperfectcomp==1) {
        if (dqdpusedeqm==0) newcarjacobians();
        else newcarjacobians_usedeqm();
    }
    
	//compute dvalue of car / de around current mpg
    for (j=1; j<=njnew; j++) dvalde1[j] = milesremaining1[j] * 0.001 * pow(mpg1[j],-2) * gas_price;
    for (j=1; j<=njnew; j++) dvalde2[j] = milesremaining2[j] * 0.001 * pow(mpg2[j],-2) * gas_price;
     
	//pull in dqdp
    for (j=1;j<=njnew;j++) {
		for (k=1;k<=njnew;k++) {
			tmpdqdp1[j][k] = dqdp1[j][k];
			tmpdqdp2[j][k] = dqdp2[j][k];
		}
	}

	//calculate dq/de
	for (j=1;j<=njnew;j++) {	
		for (k=1;k<=njnew;k++) {
			dqde1[j][k] = -dvalde1[k] * tmpdqdp1[j][k];
			dqde2[j][k] = -dvalde2[k] * tmpdqdp2[j][k];
		}
	}
	
	//first compute FOCs wrt lambdas (transformed constraints)
	//intialize
	for (m=1;m<=NMAKE;m++) {
		focs[4*njnew+m]         = 0.0;
		focs[4*njnew+NMAKE+m]   = 0.0;
		focs[4*njnew+2*NMAKE+m] = 0.0;
	}
	// and take sum of mq
	for (j=1;j<=njnew;j++) {
		focs[4*njnew+make[j]]         += mC1[j]*d1[j] + mC2[j]*d2[j];
		focs[4*njnew+NMAKE+make[j]]   += mT1[j]*d1[j] + mT2[j]*d2[j];
		focs[4*njnew+2*NMAKE+make[j]] += mP1[j]*d1[j];
	}
	
	// compute the FOCs wrt price, modular and common fuel economy
    // *** imperfect competition case ***
    if (imperfectcomp == 1) {
        for (j=1;j<=njnew;j++) {
            
            //FOCs wrt price for regions 1 and 2
            focs[j] = d1[j];
            focs[njnew+j] = d2[j];
            for (l=1;l<=4;l++) {
                lj = mc[make[j]][l]; //lj is the index for vehicle 'l' made by firm 'make'
                focs[j] += ( pparam[lj] - cost1[lj] + lambdaC[make[lj]] * mC1[lj] + lambdaT[make[lj]] * mT1[lj] + lambdaP[make[lj]] * mP1[lj] ) * tmpdqdp1[lj][j];
                focs[njnew+j] += ( pparam[njnew+lj] - cost2[lj] + lambdaC[make[lj]] * mC2[lj] + lambdaT[make[lj]] * mT2[lj] ) * tmpdqdp2[lj][j];
            }
        
                    
            //FOCs wrt fuel economy

            if (splitmpg == 0) {
                focs[2*njnew+j] = ( -dcde1[j] + lambdaC[make[j]] * dmC1[j] + lambdaT[make[j]] * dmT1[j] ) * (d1[j] + d2[j]) + lambdaP[make[j]] * dmP1[j] * d1[j];
                for (l=1;l<=4;l++) {
                    lj = mc[make[j]][l]; //lj is the index for vehicle 'l' made by firm 'make'
                    focs[2*njnew+j] += ( pparam[lj] - cost1[lj] + lambdaC[make[lj]] * mC1[lj] + lambdaT[make[lj]] * mT1[lj] + lambdaP[make[lj]] * mP1[lj] ) * dqde1[lj][j] +
                    ( pparam[njnew+lj] - cost2[lj] + lambdaC[make[lj]] * mC2[lj] + lambdaT[make[lj]] * mT2[lj]  ) * dqde2[lj][j];
                }

            }
            if (splitmpg == 1) {  //then we have two sets of focs wrt mpg
                focs[2*njnew+j] =         ( -dcde1[j] + lambdaC[make[j]]*dmC1[j] + lambdaT[make[j]]*dmT1[j] + lambdaP[make[j]]*dmP1[j] ) * d1[j] ;
                focs[4*njnew+3*NMAKE+j] = ( -dcde2[j] + lambdaC[make[j]]*dmC2[j] + lambdaT[make[j]]*dmT2[j] ) * d2[j] ;
                for (l=1;l<=4;l++) {
                    lj = mc[make[j]][l]; //lj is the index for vehicle 'l' made by firm 'make'
                    focs[2*njnew+j] += ( pparam[lj] - cost1[lj] + lambdaC[make[lj]]*mC1[lj] + lambdaT[make[lj]]*mT1[lj] + lambdaP[make[lj]]*mP1[lj] ) * dqde1[lj][j];
                    focs[4*njnew+3*NMAKE+j] += ( pparam[njnew+lj] - cost2[lj] + lambdaC[make[lj]]*mC2[lj] + lambdaT[make[lj]]*mT2[lj] ) * dqde2[lj][j];
                }
            }

            
            //FOCs wrt common fuel economy
            if (theta_1 > 0) {
                focs[3*njnew+j] = ( -1.0*dcdz1[j] * d1[j] - dcdz2[j] * d2[j] - dczdz[j]);
            }
            else {
                focs[3*njnew+j] = pparam[3*njnew+j];
            }
    
        }
    }
    
    // compute the FOCs wrt price, modular and common fuel economy
    // *** perfect competition case ***
    if (imperfectcomp == 0) {
        for (j=1;j<=njnew;j++) {
            
            //FOCs wrt price for regions 1 and 2
            //perfect competition implies the solver must choose p such that:  p - c + lambda*m = 0
  
            focs[j] = pparam[j] - cost1[j] + lambdaC[make[j]] * mC1[j] + lambdaT[make[j]] * mT1[j] + lambdaP[make[j]] * mP1[j];
            focs[njnew+j] =  pparam[njnew+j] - cost2[j] + lambdaC[make[j]] * mC2[j] + lambdaT[make[j]] * mT2[j] ;
        
                    
            //FOCs wrt fuel economy
            //perfect competition implies the solver must choose e such that:  lambda*dm/de + dval/de - dc/de = 0
            if (splitmpg == 0) {
                focs[2*njnew+j] = ( -dcde1[j] + lambdaC[make[j]] * dmC1[j] + lambdaT[make[j]] * dmT1[j] + dvalde1[j]) * (d1[j] + d2[j]) + lambdaP[make[j]] * dmP1[j] * d1[j];
            }

            if (splitmpg == 1) {  //then we have two sets of focs wrt mpg
                focs[2*njnew+j] =         ( -dcde1[j] + lambdaC[make[j]]*dmC1[j] + lambdaT[make[j]]*dmT1[j] + lambdaP[make[j]]*dmP1[j] + dvalde1[j]) * d1[j] ;
                focs[4*njnew+3*NMAKE+j] = ( -dcde2[j] + lambdaC[make[j]]*dmC2[j] + lambdaT[make[j]]*dmT2[j] + dvalde2[j]) * d2[j] ;
            }

            //FOCs wrt common fuel economy
            //common fuel economy always turned off under perfect competition (constant returns to scale)
            focs[3*njnew+j] = pparam[3*njnew+j];

        }
    }
    
	
	//scale FOCs for solver

	for (j=1;j<=njnew;j++) {
		focs[j]                 = focs[j] / 100.0;
		focs[njnew+j]           = focs[njnew+j] / 100.0;
		focs[2*njnew+j]         = focs[2*njnew+j] / 100.0;
		focs[3*njnew+j]         = focs[3*njnew+j] / 100.0;
		focs[4*njnew+3*NMAKE+j] = focs[4*njnew+3*NMAKE+j] / 100.0;
	}
	for (m=1;m<=NMAKE;m++){
		focs[4*njnew+m] = 	focs[4*njnew+m] / 100.0;
		focs[4*njnew+NMAKE+m] = focs[4*njnew+NMAKE+m] / 100.0;
		focs[4*njnew+2*NMAKE+m] = focs[4*njnew+2*NMAKE+m] / 100.0;
	}
	
	//set the FOCs for the lambda terms in cases where CAFE is not binding
	for (m=1;m<=NMAKE;m++) {
		if (cafeCstat[m] == 0) focs[4*njnew+m] = pparam[4*njnew+m];
		if (cafeTstat[m] == 0) focs[4*njnew+NMAKE+m] = pparam[4*njnew+NMAKE+m];
		if (cafePstat[m] == 0)  focs[4*njnew+2*NMAKE+m] = pparam[4*njnew+2*NMAKE+m];
	}
	
	sumfoc = 0.0;
	for (n=1;n<=numfocs;n++) sumfoc += fabs(focs[n]);
	sumfocglobal = sumfoc;


	//print for debugging
    if (verbose) {
        if (prodprob_in_jacobian == 1) printf("\nprod prob jacobian - wait...");
        else {
            printf("\n\n%lf <-ppsumfoc  foc 1 %lf foc 29 %lf foc 57 %lf lamC1 %lf p1 %lf ",sumfoc,focs[1],focs[29],focs[57],lambdaC[1],pparam[1]);
            if (dbprint) {
                for (n=1;n<=numfocs;n++)  fprintf(debug, "%20.12lf, ", pparam[n]);
                fprintf(debug,"\n");
                for (n=1;n<=numfocs;n++)  fprintf(debug, "%20.12lf, ", focs[n]);
                fprintf(debug,"\n");
                for (j=1;j<=nj;j++)  fprintf(debug, "%20.12lf, ", d1[j]);
                fprintf(debug,"\n\n");
            }
        }
    }
}


void usedcareqm(void) {
//function to solve for used vehicle rental price equilibrium given a set of new vehicle prices
	double *guess_u, *exdem_u;
	int nmkt_u; //number of markets in the used car problem
	int j, ju;
	double jacepssave;
	
	if (splitused==1) nmkt_u = njused*2; 
	else nmkt_u = njused;
	guess_u = dvector(1,nmkt_u);
	exdem_u = dvector(1,nmkt_u);
		
	//make used car rental price starting guesses (use values from last time usedcareqm was run)
	for (ju=1;ju<=njused;ju++) {
		guess_u[ju]        = rj1[njnew+ju]; //switching index from j to ju for used cars
		if (splitused == 1) {   //also make guesses for second used car market if needed
			guess_u[ju+njused] = rj2[njnew+ju]; //switching index from j to ju for used cars
		}
	}
	
	if (fullprodprob==1) {
		//call numerical solver broydn on the function returning excess demands in the used car market
		jacepssave = jaceps;
		jaceps = 0.001;
		broydn(guess_u,nmkt_u,&retvalused,usedcarexdems);
		jaceps=jacepssave;
	}
	
	//call once more to update globals (and to evaluate partial eqm demands when not in full prod prob)
	usedcarexdems(nmkt_u,guess_u,exdem_u);

	if (fullprodprob==1 & verbose) printf("\nSolved used car eqm return value %d", retvalused);
	
	free_dvector(guess_u,1,nmkt_u);	
	free_dvector(exdem_u,1,nmkt_u);
}

void usedcarexdems(int nmkt, double *depju, double *exdemu) {
	//takes in a pointer to a guess for used car prices 
	//returns excess demands in each of the nju (times 2 if split) used car markets
	
	int n, j, ju, ax, jx, jxx, tx;
	double sumed, costdiff;
    double xpscrap1[nj+1], xhj1[nj+1], xpj1[nj+1], xmilesremaining1[nj+1];
    double xpscrap2[nj+1], xhj2[nj+1], xpj2[nj+1], xmilesremaining2[nj+1];
    double fracremain1[nj+1], fracremain2[nj+1];
    double expdamagepercar[nj+1][maxt+NAGE+1], explcdamagepercar[nj+1][maxt+NAGE+1], exptau[nj+1][maxt+NAGE+1];
    double tempwavg, tempw, expfeebconstant;
           
    //trap used car rental guesses that are more negative than tau (i.e. rental plus tax is negative)
    for (ju=1;ju<=njused;ju++) {
        if (depju[ju] + (vmt[njnew+ju] / mpg1[njnew+ju]) * gas_price / 1000.0 + tau[njnew+ju]/1000.0 + proptax[njnew+ju]/1000.0 < 0.02) {
            if(verbose) printf("\nwarning: limiting negative annual cost for car %d from %lf to %lf", ju, depju[ju],0.01+0.01 * exp( depju[ju] + (vmt[njnew+ju] / mpg1[njnew+ju]) * gas_price / 1000.0 + tau[njnew+ju]/1000.0 + proptax[njnew+ju]/1000.0 - 0.02) - (vmt[njnew+ju] / mpg1[njnew+ju]) * gas_price / 1000.0 - tau[njnew+ju]/1000.0 - proptax[njnew+ju]/1000.0);
            depju[ju] = 0.01+0.01 * exp( depju[ju] + (vmt[njnew+ju] / mpg1[njnew+ju]) * gas_price / 1000.0 + tau[njnew+ju]/1000.0 + proptax[njnew+ju]/1000.0 - 0.02) - (vmt[njnew+ju] / mpg1[njnew+ju]) * gas_price / 1000.0 - tau[njnew+ju]/1000.0 - proptax[njnew+ju]/1000.0;
            boundsviolations++;
        }
    }
    
	//copy the used car prices being guessed by the solution algorithm into rj
	for (ju=1;ju<=njused;ju++) {
		rj1[njnew+ju] = depju[ju];
	}

	if (splitused == 1) {
		for (ju=1;ju<=njused;ju++) {
			rj2[njnew+ju] = depju[ju+njused];
		}
	}
	else {
		// we need to harmonize used car rental prices such that the rental plus gas price are the same when mpg differs
		for (ju=1;ju<=njused;ju++) {
			costdiff = (vmt[njnew+ju] / mpg2[njnew+ju] - vmt[njnew + ju] / mpg1[njnew+ju]) * gas_price / 1000.0;
			rj2[njnew+ju] = depju[ju] - costdiff;  //= price of region 1 car less value difference due to mpg
		}
	}

	// define expectations over future rental prices
    // if tauexpectations = mpgexpectations = 0 then expectations are fully steady-state (i.e., future depreciation expected to be the same as that implied by prices in the present time period)
    // if either of these are set larger than 0 the difference between mpg and emissions in currently-old versus future-old vehicles is used in forming expectations
    
    for (tx=t;tx<=t+18;tx++) {
        for (jx=1; jx<=nj; jx++) {
            // for time t, this loops through all cars, at t+1 loops through age 1 - 18, etc.
            expdamagepercar[jx][tx] = vmt[jx] * (damage_co * jt_phi_co[jx][tx] + damage_hc * jt_phi_hc[jx][tx] + damage_nox * jt_phi_nox[jx][tx]) / 1000000;
            explcdamagepercar[jx][tx] = vmt[jx] * (damage_co * jt_phi_co[jx][tx] + damage_hc * jt_phi_hc[jx][tx] + damage_nox * jt_phi_nox[jx][tx]) / 1000000 + proddamage[jx]*pow(0.9193,tx-t);
        }
    }

    for (tx=t;tx<=t+18;tx++) {
        
        for (jx=1+(tx-t)*njnew; jx<=nj; jx++) {
            exptau[jx][tx] = 0.0; // initialize
        }
        
        if (policytype == 1 & tx >= policystart & t >= policystart) {
            for (jx=1+(tx-t)*njnew; jx<=nj; jx++) {
                if (lcpolicy==1) exptau[jx][tx] = pollutionpolicy * explcdamagepercar[jx][tx];
                else exptau[jx][tx] = pollutionpolicy * expdamagepercar[jx][tx];
            }
        }
        
        if (policytype == 2 & tx >= policystart & t >= policystart) {
            for (ax=(tx-t); ax<=18; ax++) {
                tempwavg = 0.0;
                tempw = 0.0;
                for (jx = ax*njnew+1; jx <= (ax+1)*njnew; jx++) {
                    if (lcpolicy==1) tempwavg += explcdamagepercar[jx][tx] * based1[jx];  //current period quantity weights are expected (steady-state assumption)
                    else tempwavg += expdamagepercar[jx][tx] * based1[jx];
                    tempw += based1[jx];
                }
                tempwavg = tempwavg / tempw;
                for (jx = ax*njnew+1; jx <= (ax+1)*njnew; jx++) {
                    exptau[jx][tx] = agepolicy*tempwavg;
                }
            }
        }
        
        if (policytype == 3 & tx >= policystart & t >= policystart) {
            for (jx=1+(tx-t)*njnew; jx<=nj; jx++) {
                exptau[jx][tx] = 0.0; // note this fills all 0s in - used car tau is always 0, and new car tau doesn't need to be computed here
            }
        }
        
        if (policytype == 4 & tx >= policystart & t >= policystart) {
            for (jx=1+(tx-t)*njnew; jx<=nj; jx++) {
                exptau[jx][tx] = tau[jx]; // no change in expected tau since it is a property tax
            }
        }

        if (policytype == 5 & tx >= policystart & t >= policystart) {
            for (jx=1+(tx-t)*njnew; jx<=nj; jx++) {
                exptau[jx][tx] = tau[jx]; // no change in expected tau since it is a property tax
            }
        }
        
        if (policytype == 6 & tx >= policystart & t >= policystart) {
            tempw = 0.0;
            expfeebconstant = 0.0;
            for (jx=1; jx<=nj; jx++) {
                if (lcpolicy==1) expfeebconstant += (based1[jx]+based2[jx]) * pollutionpolicy * explcdamagepercar[jx][tx];
                else expfeebconstant += (based1[jx]+based2[jx]) * pollutionpolicy * expdamagepercar[jx][tx];
                tempw+=(based1[jx]+based2[jx]);
            }
            expfeebconstant = -expfeebconstant / tempw;

            for (jx=1+(tx-t)*njnew; jx<=nj; jx++) {
                if (lcpolicy==1) exptau[jx][tx] = expfeebconstant + pollutionpolicy * explcdamagepercar[jx][tx];
                else exptau[jx][tx] = expfeebconstant + pollutionpolicy * expdamagepercar[jx][tx];
            }
        }
        
        if (policytype == 7 & tx >= policystart & t >= policystart) {
            tempw = 0.0;
            expfeebconstant = 0.0;
            for (jx=1; jx<=nj; jx++) {
                if (age[jx] >= policyagestart) expfeebconstant += (based1[jx]+based2[jx]) * pollutionpolicy;
                tempw+=(based1[jx]+based2[jx]);
            }
            expfeebconstant = -expfeebconstant / tempw;

            for (jx=1+(tx-t)*njnew; jx<=nj; jx++) {
                if (age[jx] >= policyagestart) exptau[jx][tx] = expfeebconstant + pollutionpolicy;
                else exptau[jx][tx] = 0.0;
            }
        }
        
        
        if (policytype == 8 & tx >= policystart & t >= policystart) {

            for (jx=1+(tx-t)*njnew; jx<=nj; jx++) {
                if (age[jx] >= policyagestart) exptau[jx][tx] = pollutionpolicy;
                else exptau[jx][tx] = 0.0;
            }
        }
        

        if (policytype == 10 & tx >= policystart & t >= policystart) {
            for (jx=1+(tx-t)*njnew; jx<=nj; jx++) {
                exptau[jx][tx] = tau_read[jx][tx];
            }
        }
    }
    
    //now can adjust expected rental prices
    for (tx=t;tx<=t+18;tx++) {
        for (jx=1+(tx-t)*njnew; jx<=nj; jx++) {
            //benchmark expectation is current rental price of same-age vehicle, i.e. steady-state
            exprj1[jx][tx] = rj1[jx];
            exprj2[jx][tx] = rj2[jx];
   
            //increase the expected rental price to account for gasoline savings relative to current versions of the same vehicle:
            if (mpgexpectations>0) {
                exprj1[jx][tx] += mpgexpectations * 0.001*vmt[jx]*gas_price * (1/mpg1[jx] - 1/mpg1[jx-(tx-t)*njnew]);
                exprj2[jx][tx] += mpgexpectations * 0.001*vmt[jx]*gas_price * (1/mpg2[jx] - 1/mpg2[jx-(tx-t)*njnew]);
            }
            
            //increase the expected rental price to account for emissions tax savings relative to current versions of the same vehicle:
            if (tauexpectations>0) {
                exprj1[jx][tx] +=  tauexpectations * (exptau[jx][t] - exptau[jx][tx]) / 1000;
                exprj2[jx][tx] +=  tauexpectations * (exptau[jx][t] - exptau[jx][tx]) / 1000;
            }
        }
    }

    // determine the final price of vehicle j
    for (j=1;j<=nj;j++) {
        
        //loop through from future oldest "jx" versions of car j backward
        for (ax=18; ax>=age[j]; ax--) {
            jx = j+(ax-age[j])*njnew;
            if (ax==18) {
                xpj1[jx] = exprj1[jx][t+ax-age[j]]; //this is the price we expect car j to have when it is 18 years old (and so has become car "jx")
                xpj2[jx] = exprj2[jx][t+ax-age[j]];
                xmilesremaining1[jx] = vmt[jx];
                xmilesremaining2[jx] = vmt[jx];
            }
            if (ax>0 && ax<18) {
                xpj1[jx] = exprj1[jx][t+ax-age[j]] + (1.0 - xpscrap1[jx+njnew])*(xpj1[jx+njnew]-xhj1[jx+njnew])/(1.0 + drate);
                xpj2[jx] = exprj2[jx][t+ax-age[j]] + (1.0 - xpscrap2[jx+njnew])*(xpj2[jx+njnew]-xhj2[jx+njnew])/(1.0 + drate);
                xmilesremaining1[jx] = vmt[jx] + (1.0 - xpscrap1[jx+njnew]) * xmilesremaining1[jx+njnew] / (1.0 + drate);
                xmilesremaining2[jx] = vmt[jx] + (1.0 - xpscrap2[jx+njnew]) * xmilesremaining2[jx+njnew] / (1.0 + drate);
            }
            
            //new cars
            if (ax==0) {
                //determine rj for new cars; note pj for new cars is already populated
                if (rj_expectations_static==0) {  // then get implicit rj for new cars
                    exprj1[jx][t] = pj1[jx] - (1.0 - xpscrap1[jx+njnew]) * (xpj1[jx+njnew]-xhj1[jx+njnew])/(1.0 + drate);
                    exprj2[jx][t] = pj2[jx] - (1.0 - xpscrap2[jx+njnew]) * (xpj2[jx+njnew]-xhj2[jx+njnew])/(1.0 + drate);
                }
                if (rj_expectations_static==1) {  // then use a fixed ratio of first-period rental to purchase price
                    exprj1[jx][t] = pj1[jx] * baserj1[jx] / basepj1[jx];
                    exprj2[jx][t] = pj2[jx] * baserj2[jx] / basepj2[jx];
                }
                xmilesremaining1[jx] = vmt[jx] + (1.0 - xpscrap1[jx+njnew]) * xmilesremaining1[jx+njnew] / (1.0 + drate);
                xmilesremaining2[jx] = vmt[jx] + (1.0 - xpscrap2[jx+njnew]) * xmilesremaining2[jx+njnew] / (1.0 + drate);
                
       
                //trap implied new car rentals that are more negative than tau and gas expend (i.e. negative ownership cost)
                if (exprj1[jx][t] + (vmt[jx] / mpg1[jx]) * gas_price / 1000.0 + tau[jx]/1000.0 + proptax[jx]/1000.0 < 0.02) {
                    if(verbose) printf("\nwarning: negative implied rental for new car %d", jx);
                    exprj1[jx][t] = 0.01+0.01 * exp( exprj1[jx][t] + (vmt[jx] / mpg1[jx]) * gas_price / 1000.0 + tau[jx]/1000.0 + proptax[jx]/1000.0 - 0.02) - (vmt[jx] / mpg1[jx]) * gas_price / 1000.0 - tau[jx]/1000.0 - proptax[jx]/1000.0;
                    exprj2[jx][t] = exprj1[jx][t];
                    boundsviolations++;
                }
            
                //value of new cars to the rental company (calculated as above, this will equal price if rj_expectations_static=0):
                valj1[jx] = exprj1[jx][t] + (1.0 - xpscrap1[jx+njnew])* (xpj1[jx+njnew]-xhj1[jx+njnew])/(1.0 + drate);
                valj2[jx] = exprj2[jx][t] + (1.0 - xpscrap2[jx+njnew])* (xpj2[jx+njnew]-xhj2[jx+njnew])/(1.0 + drate);

                xpscrap1[jx] = 0.0;
                xhj1[jx] = 0.0;
                xpscrap2[jx] = 0.0;
                xhj2[jx] = 0.0;
                
                //payment model
                fracremain1[j] = 1.0;  //fraction of each new car expected to remain over time (age in this notation)
                fracremain2[j] = 1.0;
                for (jxx=j+njnew; jxx<=nj; jxx+=njnew) {
                    fracremain1[jxx] = fracremain1[jxx-njnew] * (1-xpscrap1[jxx]);
                    fracremain2[jxx] = fracremain2[jxx-njnew] * (1-xpscrap2[jxx]);
                }
                for (jxx=j; jxx<=nj; jxx+=njnew) {
                    //payments on 1 year old cars for next year (and 2 year old cars in 2 years, etc.)
                    payment1[jxx][t+age[jxx]] = fracremain1[jxx] * ( exprj1[jxx][t+age[jxx]] - xhj1[jxx]);
                    payment2[jxx][t+age[jxx]] = fracremain2[jxx] * ( exprj2[jxx][t+age[jxx]] - xhj2[jxx]);

                }
            }
            
            //determine scrap rate
            if (ax>0) {
                if (xpj1[jx]>0) {
                    xpscrap1[jx] = b[jx] * pow(xpj1[jx],eta[jx]);
                    xhj1[jx] = (-1.0 * pow(b[jx],-1/eta[jx]) * eta[jx] + b[jx] * eta[jx] * pow(xpj1[jx],1+eta[jx])) /
                    ((1+eta[jx]) * (-1.0 + b[jx] * pow(xpj1[jx],eta[jx])));
                }
                else {
                    xpscrap1[jx] = 1.0;
                    xhj1[jx] = 0.0;
                    boundsviolations++;
                }
                if (xpj2[jx]>0) {
                    xpscrap2[jx] = b[jx] * pow(xpj2[jx],eta[jx]);
                    xhj2[jx] = (-1.0 * pow(b[jx],-1/eta[jx]) * eta[jx] + b[jx] * eta[jx] * pow(xpj2[jx],1+eta[jx])) /
                    ((1+eta[jx]) * (-1.0 + b[jx] * pow(xpj2[jx],eta[jx])));
                }
                else {
                    xpscrap2[jx] = 1.0;
                    xhj2[jx] = 0.0;
                    boundsviolations++;
                }
                
                xpscrap1[jx] = DMIN(xpscrap1[jx],0.9999);
                xpscrap2[jx] = DMIN(xpscrap2[jx],0.9999);
            }
            
            //if jx=j then we've finished working our way backward through future ages and can assign values to the current-period arrays
            if (jx==j) {
                if (ax>0) pj1[j] = xpj1[j];
                else rj1[j] = exprj1[j][t];
                pscrap1[j] = xpscrap1[j];
                hj1[j] = xhj1[j];
                milesremaining1[j] = xmilesremaining1[j];
                if (ax>0) pj2[j] = xpj2[j];
                else rj2[j] = exprj2[j][t];
                pscrap2[j] = xpscrap2[j];
                hj2[j] = xhj2[j];
                milesremaining2[j] = xmilesremaining2[j];
            }
            if (jx==j && ax<18) {
                //save current expectations about prices etc. one year forward (for use in capital gains calculations)
                exppj1[j+njnew] = xpj1[j+njnew];
                exppscrap1[j+njnew] = xpscrap1[j+njnew];
                exphj1[j+njnew] = xhj1[j+njnew];
                exppj2[j+njnew] = xpj2[j+njnew];
                exppscrap2[j+njnew] = xpscrap2[j+njnew];
                exphj2[j+njnew] = xhj2[j+njnew];
            }
        }
    }
        
    //run the household problem:
    hhproblem(rj1, rj2, d1, d2, M_1, M_2);

    // Run it a few more times in policy cases so the income adjustments from profits (which are a function of demand) converge.  Note this won't affect the convergence criteria or solution, but it can speed the algorithm since it will improve accuracy of the numerical derivatives.
    if (runmode==1) {
        for(n=1;n<=3;n++) hhproblem(rj1, rj2, d1, d2, M_1, M_2);
    }
	
	//excess demands in the njused*2 markets are:
	for (ju=1;ju<=njused;ju++) {
		exdemu[ju] = d1[njnew+ju] - (1.0 - pscrap1[njnew+ju]) * maxsup1[njnew+ju];
		if (splitused == 1) {
			exdemu[ju+njused] = d2[njnew+ju] - (1.0 - pscrap2[njnew+ju]) * maxsup2[njnew+ju];
		}
		else {
			exdemu[ju] += d2[njnew+ju] - (1.0 - pscrap2[njnew+ju]) * maxsup2[njnew+ju];
		}

		//scale
		exdemu[ju] = exdemu[ju]/100.0;
		if (splitused == 1) exdemu[ju+njused] = exdemu[ju+njused]/1000.0;
	}

	//print some values for debugging
	if (usedcar_in_jacobian == 0) {
		sumed = 0.0;
		for (ju=1;ju<=nmkt;ju++)  sumed += fabs(exdemu[ju]);
		if (verbose) printf(" SUMED %lf", sumed);
		sumedglobal = sumed;
		
		//prints all guesses and ed's to the debug file
        if (dbprint) {
            for (j=1;j<=njnew;j++)  fprintf(debug, "%20.12lf, ", rj1[j]);
            for (n=1;n<=nmkt;n++)  fprintf(debug, "%20.12lf, ", depju[n]);
            fprintf(debug,"\n");
            for (j=1;j<=nj;j++)  fprintf(debug, "%20.12lf, ", pj1[j]);
            fprintf(debug,"\n");
            for (j=1;j<=nj;j++)  fprintf(debug, "%20.12lf, ", d1[j] + d2[j]);
            fprintf(debug,"\n");
             for (j=1;j<=nj;j++)  fprintf(debug, "%20.12lf, ", pscrap1[j]);
            fprintf(debug,"\n");
            for (j=1;j<=njnew;j++) fprintf(debug, " , ");
            for (n=1;n<=nmkt;n++)  fprintf(debug, "%20.12lf, ", exdemu[n]);
            fprintf(debug,"\n\n");
        }
	}
}

void hhproblem(double *hhrj1, double *hhrj2, double *hhd1, double *hhd2, double hhM_1, double hhM_2) {
	//solves the households' problems. takes rental prices as an input and returns demands as outputs
    //makes use of globals (e.g. pj and valj) in updating profits for policy runs only
	//the calibrated alpha values are outputs from the calibration program
	//all composite prices are normalized to 1.
	//using this data, the subroutine first computes demand ratios at the lowest nest (nest 1), moving upwards to nest 4
	//the subroutine then computes demand levels at the highest nest (nest 5)
	//finally, working backwards, the subroutine computes the demand levels for lower nests
	int j;
	int dim1;
	int dim2;
	int dim3;
	int dim4;
	double temp1;
	double temp2;
    double tempw;
      
    //INCOME ADJUSTMENTS
    //note that these rely on prior evaluation demands (since income and demand are interdependent)
    //there are also no income adjustments in the base case (profits are assumed to be part of exogenously read gdp)
    
    //compute new car producer profits (using demands from previous hhprob evaluation, prices and costs are from current guess in prodprob)
    profit1 = 0.0;
    profit2 = 0.0;
    
    for (j=1;j<=njnew;j++) {
        profit1 += (pj1[j] - cost1[j]) * d1[j] - (M_total_1 / (M_total_1 + M_total_2)) * costz[j];
        profit2 += (pj2[j] - cost2[j]) * d2[j] - (M_total_2 / (M_total_1 + M_total_2)) * costz[j];
    }
    

    if (depreciation_mode == 1) {
        //full depreciation of new car purchases and repairs in current time period
        //profit of rental company is: "current rental income, less the cost of repairs and replacements to the fleet"
        profitrental1 = 0.0;
        profitrental2 = 0.0;
        for (j=1;j<=nj;j++) {   //rental income less current repair spending
            profitrental1 += d1[j] * (rj1[j] - hj1[j]);
            profitrental2 += d2[j] * (rj2[j] - hj2[j]);
        }
        for (j=1;j<=njnew;j++) {   //less cost of replacement vehicles
            profitrental1 -= d1[j] * pj1[j];
            profitrental2 -= d2[j] * pj2[j];
        }
    }
    
    
    if (depreciation_mode == 2) {
        //excess depreciation (i.e. capital loss) between periods is charged against profits
        //profit of rental company is: "current rental income, less current expected depreciation, less capital loss between time periods
        profitrental1 = 0.0;
        profitrental2 = 0.0;
        for (j=njnew+1;j<=nj;j++) {
            profitrental1 += lastd1[j-njnew]*( (1-pscrap1[j])*(pj1[j]-hj1[j]) - (1-lastexppscrap1[j])*(lastexppj1[j]-lastexphj1[j]) );
            profitrental2 += lastd2[j-njnew]*( (1-pscrap2[j])*(pj2[j]-hj2[j]) - (1-lastexppscrap2[j])*(lastexppj2[j]-lastexphj2[j]) );
        }
    }
    
    
    if (depreciation_mode == 3) {
        //profit is current rentals less repair expenditure, less fixed payments made on past vehicle purchases
        //note that for age 0 vehicles the rental and payment are equal, and there are no repairs, so they don't appear in profits
        profitrental1 = 0.0;
        profitrental2 = 0.0;
        for (j=njnew+1;j<=nj;j++) {   //rental income less current repair spending
            profitrental1 += d1[j] * (rj1[j] - hj1[j]);
            profitrental2 += d2[j] * (rj2[j] - hj2[j]);
        }
        for (j=njnew+1; j<=nj; j++) {  //payments
            profitrental1 -= d1_when_new[j] * payment1[j][t];
            profitrental2 -= d2_when_new[j] * payment2[j][t];
        }
    }
    
    //update tau when it depends on pj (no need to update for most cases, where it depends on e.g. phi)
    if (policytype == 5 & t >= policystart) {
        for (j=1; j<=nj; j++) {
            tau[j] = 2 * proptaxpolicy * pj1[j] * 1000;
        }
    }
    
    //calculate total tax bill (from the policy taxes, tau) to add as lump-sum transfer
    taxrev_tot_1 = 0.0;
    taxrev_tot_2 = 0.0;
    for (j = 1; j <= nj; j++) {
        taxrev_tot_1 += tau[j] * d1[j] / 1000.0;
        taxrev_tot_2 += tau[j] * d2[j] / 1000.0;
    }

    
    //calculate the property tax revenue to add (relative to baseline) as a transfer:
    proptaxrev1 = 0.0;
    proptaxrev2 = 0.0;
    for (j = 1; j <= nj; j++) {
        proptax[j] = 2 * proptaxrate * pj1[j] * 1000;
        proptaxrev1 += proptax[j] * d1[j] / 1000.0;
        proptaxrev2 += proptax[j] * d2[j] / 1000.0;
    }
    
    if (runmode==0) {
        baseprofit1 = profit1;
        baseprofit2 = profit2;
        baseprofitrental1 = profitrental1;
        baseprofitrental2 = profitrental2;
        baseproptaxrev1 = proptaxrev1;
        baseproptaxrev2 = proptaxrev2;
    }

    //share out the change in profits (new and used) by region, and add to total income
    profit_adj1 = (profit1 + profit2 - baseprofit1 - baseprofit2 + profitrental1 + profitrental2 - baseprofitrental1 - baseprofitrental2) * (hhM_1 / (hhM_1 + hhM_2));
	profit_adj2 = (profit1 + profit2 - baseprofit1 - baseprofit2 + profitrental1 + profitrental2 - baseprofitrental1 - baseprofitrental2) * (hhM_2 / (hhM_1 + hhM_2));
	

    M_total_1 = hhM_1 + profit_adj1 + taxrev_tot_1 + proptaxrev1 - baseproptaxrev1;
    M_total_2 = hhM_2 + profit_adj2 + taxrev_tot_2 + proptaxrev2 - baseproptaxrev2;
    
    if (policytype==4 & t >= policystart) {   // revenue neutral policy for tau (i.e. "taxrev_tot" variable will be zero in eqm, so don't add in out of eqm values)
        M_total_1 = hhM_1 + profit_adj1 + proptaxrev1 - baseproptaxrev1;
        M_total_2 = hhM_2 + profit_adj2 + proptaxrev2 - baseproptaxrev2;
        
        if(d1[1] > 0) {
        tempw = 0.0;
        proptaxconstant = 0.0;  // update constant with latest demands
        for (j=1; j<=nj; j++) {
            proptaxconstant += (d1[j]+d2[j]) * 2* proptaxpolicy * pj1[j] * 1000;
            tempw+=(d1[j]+d2[j]);
        }
        proptaxconstant = -proptaxconstant / tempw;
        for (j=1; j<=nj; j++) {
            tau[j] = proptaxconstant + 2* proptaxpolicy * pj1[j] * 1000;
        }
        }
         
    }
    
    if (policytype==6 & t >= policystart) {   // revenue neutral policy for tau (i.e. "taxrev_tot" variable will be zero in eqm, so don't add in out of eqm values)
        M_total_1 = hhM_1 + profit_adj1 + proptaxrev1 - baseproptaxrev1;
        M_total_2 = hhM_2 + profit_adj2 + proptaxrev2 - baseproptaxrev2;
        
        if(d1[1] > 0) {
            tempw = 0.0;
            feebconstant = 0.0;  // update constant with latest demands
            for (j=1; j<=nj; j++) {
                if (lcpolicy==1) feebconstant += (d1[j]+d2[j]) * pollutionpolicy * lcdamagepercar[j];
                else feebconstant += (d1[j]+d2[j]) * pollutionpolicy * damagepercar[j];
                tempw+=(d1[j]+d2[j]);
            }
            feebconstant = -feebconstant / tempw;
            for (j=1; j<=nj; j++) {
                if (lcpolicy==1) tau[j] = feebconstant +  pollutionpolicy * lcdamagepercar[j];
                else tau[j] = feebconstant +  pollutionpolicy * damagepercar[j];
            }
        }
    }
    
    
    if (policytype==7 & t >= policystart) {   // revenue neutral policy for tau (i.e. "taxrev_tot" variable will be zero in eqm, so don't add in out of eqm values)
        M_total_1 = hhM_1 + profit_adj1 + proptaxrev1 - baseproptaxrev1;
        M_total_2 = hhM_2 + profit_adj2 + proptaxrev2 - baseproptaxrev2;
        
        if(d1[1] > 0) {
            tempw = 0.0;
            feebconstant = 0.0;  // update constant with latest demands
            for (j=1; j<=nj; j++) {
                if (age[j]>=policyagestart) feebconstant += (d1[j]+d2[j]) * pollutionpolicy;
                tempw+=(d1[j]+d2[j]);
            }
            feebconstant = -feebconstant / tempw;
            for (j=1; j<=nj; j++) {
                if (age[j]>=policyagestart) tau[j] = feebconstant +  pollutionpolicy;
                else tau[j]=0.0;
            }
        }
    }
    
    
    if (M_total_1 < 1) {
        M_total_1 = 1;
        boundsviolations++;
    }
    if (M_total_2 < 1) {
        M_total_2 = 1;
        boundsviolations++;
    }
    
    //DEMAND FUNCTION:
	//copy new values for rental prices (in hhrj1 and hhrj2) into the tap_1_1 and tap_1_2 arrays
	//and add the fuel cost.
	for (dim1=1; dim1<=19; dim1++) {
		for (dim2=1; dim2<=2; dim2++) {
			for (dim3=1; dim3<=2; dim3++) {
				for (dim4=1; dim4<=7; dim4++) {
					tap_1_1[dim1][dim2][dim3][dim4] =
                    hhrj1[((dim1 - 1)*(NTYPE * NSIZE * NMAKE) + (dim2 - 1)*(NSIZE * NMAKE) + (dim3 - 1)*NMAKE + dim4)]
                    + (vmt[((dim1 - 1)*(NTYPE * NSIZE * NMAKE) + (dim2 - 1)*(NSIZE * NMAKE) + (dim3 - 1)*NMAKE + dim4)] / mpg1[((dim1 - 1)*(NTYPE * NSIZE * NMAKE) + (dim2 - 1)*(NSIZE * NMAKE) + (dim3 - 1)*NMAKE + dim4)]) * gas_price / 1000.0
                    + tau[((dim1 - 1)*(NTYPE * NSIZE * NMAKE) + (dim2 - 1)*(NSIZE * NMAKE) + (dim3 - 1)*NMAKE + dim4)] / 1000.0
                    + proptax[((dim1 - 1)*(NTYPE * NSIZE * NMAKE) + (dim2 - 1)*(NSIZE * NMAKE) + (dim3 - 1)*NMAKE + dim4)] / 1000.0;
					tap_1_2[dim1][dim2][dim3][dim4] =
                    hhrj2[((dim1 - 1)*(NTYPE * NSIZE * NMAKE) + (dim2 - 1)*(NSIZE * NMAKE) + (dim3 - 1)*NMAKE + dim4)]
                    + (vmt[((dim1 - 1)*(NTYPE * NSIZE * NMAKE) + (dim2 - 1)*(NSIZE * NMAKE) + (dim3 - 1)*NMAKE + dim4)] / mpg2[((dim1 - 1)*(NTYPE * NSIZE * NMAKE) + (dim2 - 1)*(NSIZE * NMAKE) + (dim3 - 1)*NMAKE + dim4)]) * gas_price / 1000.0
                    + tau[((dim1 - 1)*(NTYPE * NSIZE * NMAKE) + (dim2 - 1)*(NSIZE * NMAKE) + (dim3 - 1)*NMAKE + dim4)] / 1000.0
                    + proptax[((dim1 - 1)*(NTYPE * NSIZE * NMAKE) + (dim2 - 1)*(NSIZE * NMAKE) + (dim3 - 1)*NMAKE + dim4)] / 1000.0;
                    if (tap_1_1[dim1][dim2][dim3][dim4] < 0.0001) tap_1_1[dim1][dim2][dim3][dim4] = 0.0001;
                    if (tap_1_2[dim1][dim2][dim3][dim4] < 0.0001) tap_1_2[dim1][dim2][dim3][dim4] = 0.0001;
				}
			}
		}
	}
	
	// Initialize all (nest 1) demands to 0
	for (j=0;j<=nj;j++) {
		hhd1[j] = 0.0;
		hhd2[j] = 0.0;
	}
	
	// Compute composite prices for nest 1
	for (dim2=1; dim2<=2; dim2++) {
		for (dim3=1; dim3<=2; dim3++) {
			for (dim1=1; dim1<=19; dim1++) {
				temp1 = 0.0; 
				temp2 = 0.0;
				for (dim4=1; dim4<=7; dim4++) {
					temp1 += pow(alpha_1_1[dim1][dim2][dim3][dim4], (1.0 / (1.0 - rho_1[dim1][dim2][dim3]))) 
					* pow(tap_1_1[dim1][dim2][dim3][dim4], (rho_1[dim1][dim2][dim3] / (rho_1[dim1][dim2][dim3] - 1.0)));
					temp2 += pow(alpha_1_2[dim1][dim2][dim3][dim4], (1.0 / (1.0 - rho_1[dim1][dim2][dim3]))) 
					* pow(tap_1_2[dim1][dim2][dim3][dim4], (rho_1[dim1][dim2][dim3] / (rho_1[dim1][dim2][dim3] - 1.0)));
				}
				tap_comp_1_1[dim1][dim2][dim3] = pow(temp1, ((rho_1[dim1][dim2][dim3] - 1.0) / rho_1[dim1][dim2][dim3]));
				tap_comp_1_2[dim1][dim2][dim3] = pow(temp2, ((rho_1[dim1][dim2][dim3] - 1.0) / rho_1[dim1][dim2][dim3]));
			}
		}
	}
	
	//compute demand ratios for nest 1
	for (dim2=1; dim2<=2; dim2++) {
		for (dim3=1; dim3<=2; dim3++) {
			for (dim1=1; dim1<=19; dim1++) {
				for (dim4=1; dim4<=7; dim4++) {
					ratio_1_1[dim1][dim2][dim3][dim4] = pow(((alpha_1_1[dim1][dim2][dim3][dim4] * tap_comp_1_1[dim1][dim2][dim3]) / tap_1_1[dim1][dim2][dim3][dim4]), (1 / (1 - rho_1[dim1][dim2][dim3])));
					ratio_1_2[dim1][dim2][dim3][dim4] = pow(((alpha_1_2[dim1][dim2][dim3][dim4] * tap_comp_1_2[dim1][dim2][dim3]) / tap_1_2[dim1][dim2][dim3][dim4]), (1 / (1 - rho_1[dim1][dim2][dim3])));
				}
			}
		}
	}
	
	//compute composite prices for nest 2
	for (dim2=1; dim2<=2; dim2++) {
		for (dim3=1; dim3<=2; dim3++) {
			temp1 = 0.0; 
			temp2 = 0.0;
			for (dim1=1; dim1<=19; dim1++) {
				temp1 += pow(alpha_2_1[dim1][dim2][dim3], (1.0 / (1.0 - rho_2[dim2][dim3]))) 
				* pow(tap_comp_1_1[dim1][dim2][dim3], (rho_2[dim2][dim3] / (rho_2[dim2][dim3] - 1.0)));
				temp2 += pow(alpha_2_2[dim1][dim2][dim3], (1.0 / (1.0 - rho_2[dim2][dim3]))) 
				* pow(tap_comp_1_2[dim1][dim2][dim3], (rho_2[dim2][dim3] / (rho_2[dim2][dim3] - 1.0)));
			}
			tap_comp_2_1[dim2][dim3] = pow(temp1, ((rho_2[dim2][dim3] - 1.0) / rho_2[dim2][dim3]));
			tap_comp_2_2[dim2][dim3] = pow(temp2, ((rho_2[dim2][dim3] - 1.0) / rho_2[dim2][dim3]));
		}
	}
	
	//compute demand ratios for nest 2
	for (dim2=1; dim2<=2; dim2++) {
		for (dim3=1; dim3<=2; dim3++) {
			for (dim1=1; dim1<=19; dim1++) {
				ratio_2_1[dim1][dim2][dim3] = pow(((alpha_2_1[dim1][dim2][dim3] * tap_comp_2_1[dim2][dim3]) / tap_comp_1_1[dim1][dim2][dim3]), (1 / (1 - rho_2[dim2][dim3])));
				ratio_2_2[dim1][dim2][dim3] = pow(((alpha_2_2[dim1][dim2][dim3] * tap_comp_2_2[dim2][dim3]) / tap_comp_1_2[dim1][dim2][dim3]), (1 / (1 - rho_2[dim2][dim3])));
			}
		}
	}
	
	//compute composite prices for nest 3
	for (dim2=1; dim2<=2; dim2++) {
		temp1 = 0.0; 
		temp2 = 0.0;
		for (dim3=1; dim3<=2; dim3++) {
			temp1 += pow(alpha_3_1[dim2][dim3], (1.0 / (1.0 - rho_3[dim2]))) 
			* pow(tap_comp_2_1[dim2][dim3], (rho_3[dim2] / (rho_3[dim2] - 1.0)));
			temp2 += pow(alpha_3_2[dim2][dim3], (1.0 / (1.0 - rho_3[dim2]))) 
			* pow(tap_comp_2_2[dim2][dim3], (rho_3[dim2] / (rho_3[dim2] - 1.0)));
		}
		tap_comp_3_1[dim2] = pow(temp1, ((rho_3[dim2] - 1.0) / rho_3[dim2]));
		tap_comp_3_2[dim2] = pow(temp2, ((rho_3[dim2] - 1.0) / rho_3[dim2]));
	}
	
	//compute demand ratios for nest 3
	for (dim2=1; dim2<=2; dim2++) {
		for (dim3=1; dim3<=2; dim3++) {
			ratio_3_1[dim2][dim3] = pow(((alpha_3_1[dim2][dim3] * tap_comp_3_1[dim2]) / tap_comp_2_1[dim2][dim3]), (1 / (1 - rho_3[dim2])));
			ratio_3_2[dim2][dim3] = pow(((alpha_3_2[dim2][dim3] * tap_comp_3_2[dim2]) / tap_comp_2_2[dim2][dim3]), (1 / (1 - rho_3[dim2])));
		}
	}
	
	//compute composite prices for nest 4
	temp1 = 0.0; 
	temp2 = 0.0;
	for (dim2=1; dim2<=2; dim2++) {
		temp1 += pow(alpha_4_1[dim2], (1.0 / (1.0 - rho_4))) 
		* pow(tap_comp_3_1[dim2], (rho_4 / (rho_4 - 1.0)));
		temp2 += pow(alpha_4_2[dim2], (1.0 / (1.0 - rho_4))) 
		* pow(tap_comp_3_2[dim2], (rho_4 / (rho_4 - 1.0)));
	}
	tap_comp_4_1 = pow(temp1, ((rho_4 - 1.0) / rho_4));
	tap_comp_4_2 = pow(temp2, ((rho_4 - 1.0) / rho_4));
	
	//compute demand ratios for nest 4
	for (dim2=1; dim2<=2; dim2++) {
		ratio_4_1[dim2] = pow(((alpha_4_1[dim2] * tap_comp_4_1) / tap_comp_3_1[dim2]), (1 / (1 - rho_4)));
		ratio_4_2[dim2] = pow(((alpha_4_2[dim2] * tap_comp_4_2) / tap_comp_3_2[dim2]), (1 / (1 - rho_4)));
	}
	
	//compute demand levels for nest 5
	temp1 = pow(alpha_5_1[1],(1 /(1 - rho_5))) * pow(tap_comp_4_1,(rho_5 /(rho_5 - 1))) + pow(alpha_5_1[2],(1 /(1 - rho_5))) * pow(tap_x_1,(rho_5 /(rho_5 - 1)));
	temp2 = pow(alpha_5_2[1],(1 /(1 - rho_5))) * pow(tap_comp_4_2,(rho_5 /(rho_5 - 1))) + pow(alpha_5_2[2],(1 /(1 - rho_5))) * pow(tap_x_2,(rho_5 /(rho_5 - 1)));
	//demand for composite vehicles in region 1
	dem_comp_5_1[1] = pow((alpha_5_1[1] / tap_comp_4_1),(1 /(1 - rho_5))) * (M_total_1 / temp1);
	//demand for other goods in region 1
	dem_comp_5_1[2] = pow((alpha_5_1[2] / tap_x_1),(1 /(1 - rho_5))) * (M_total_1 / temp1);
	//demand for composite vehicles in region 2
	dem_comp_5_2[1] = pow((alpha_5_2[1] / tap_comp_4_2),(1 /(1 - rho_5))) * (M_total_2 / temp2);
	//demand for other goods in region 2
	dem_comp_5_2[2] = pow((alpha_5_2[2] / tap_x_2),(1 /(1 - rho_5))) * (M_total_2 / temp2);
	
	//compute demand levels for nest 4
	for (dim2=1; dim2<=2; dim2++) {
		dem_comp_4_1[dim2] = ratio_4_1[dim2] * dem_comp_5_1[1];
		dem_comp_4_2[dim2] = ratio_4_2[dim2] * dem_comp_5_2[1];
	}
	
	//compute demand levels for nest 3
	for (dim2=1; dim2<=2; dim2++) {
		for (dim3=1; dim3<=2; dim3++) {
			dem_comp_3_1[dim2][dim3] = ratio_3_1[dim2][dim3] * dem_comp_4_1[dim2];
			dem_comp_3_2[dim2][dim3] = ratio_3_2[dim2][dim3] * dem_comp_4_2[dim2];
		}
	}
	
	//compute demand levels for nest 2
	for (dim2=1; dim2<=2; dim2++) {
		for (dim3=1; dim3<=2; dim3++) {
			for (dim1=1; dim1<=19; dim1++) {
				dem_comp_2_1[dim1][dim2][dim3] = ratio_2_1[dim1][dim2][dim3] * dem_comp_3_1[dim2][dim3];
				dem_comp_2_2[dim1][dim2][dim3] = ratio_2_2[dim1][dim2][dim3] * dem_comp_3_2[dim2][dim3];
			}
		}
	}
	
	//compute demand levels for nest 1
    //if (usedcar_in_jacobian == 0) fprintf(debug2,"\n");
	for (dim2=1; dim2<=2; dim2++) {
		for (dim3=1; dim3<=2; dim3++) {
			for (dim1=1; dim1<=19; dim1++) {
				for (dim4=1; dim4<=7; dim4++) {
					dem_1_1[dim1][dim2][dim3][dim4] = ratio_1_1[dim1][dim2][dim3][dim4] * dem_comp_2_1[dim1][dim2][dim3];
					dem_1_2[dim1][dim2][dim3][dim4] = ratio_1_2[dim1][dim2][dim3][dim4] * dem_comp_2_2[dim1][dim2][dim3];
                    //if (usedcar_in_jacobian == 0) fprintf(debug2,"%lf,",dem_1_1[dim1][dim2][dim3][dim4]);
				}
			}
		}
	}
		
	
	for (dim1=1; dim1<=19; dim1++) {
		for (dim2=1; dim2<=2; dim2++) {
			for (dim3=1; dim3<=2; dim3++) {
				for (dim4=1; dim4<=7; dim4++) {
					hhd1[((dim1 - 1)*(NTYPE * NSIZE * NMAKE) + (dim2 - 1)*(NSIZE * NMAKE) + (dim3 - 1)*NMAKE + dim4)] = dem_1_1[dim1][dim2][dim3][dim4];
					hhd2[((dim1 - 1)*(NTYPE * NSIZE * NMAKE) + (dim2 - 1)*(NSIZE * NMAKE) + (dim3 - 1)*NMAKE + dim4)] = dem_1_2[dim1][dim2][dim3][dim4];
				}
			}
		}
	}
	
}

void indirectutil(double v_1, double x_1, double v_2, double x_2) {
	util1[t][runmode] = pow((alpha_5_1[1] * pow(v_1,rho_5) + alpha_5_1[2] * pow(x_1,rho_5)) , (1.0 / rho_5));
	util2[t][runmode] = pow((alpha_5_2[1] * pow(v_2,rho_5) + alpha_5_2[2] * pow(x_2,rho_5)) , (1.0 / rho_5));
}


void assign_bindingness(int t)	{
//sets the fuel economy limits for time period t
	int m;

	for (m=1; m<=NMAKE; m++) {

        cafelimC[m] = cafeCtarget[t];
        cafelimT[m] = cafeTtarget[t];

        // use the calibrated lambdas in the base version of time 1, combined with the time 1 base CAFE values, to adjust the limits if there is inconsistency:
        // if binding in time 1 (as indicated by lambda), and there is a discrepancy between the time 1 limit as read in and the actual time 1 value, then shift the limit to reconcile:
        if (lambdaCread[m]>0) cafelimC[m] = cafeCtarget[t] + cafeC_baset1[m]-cafeCtarget[1];
        if (lambdaTread[m]>0) cafelimT[m] = cafeTtarget[t] + cafeT_baset1[m]-cafeTtarget[1];
        // if not binding in time 1 (as indicated by lambda), but the time 1 limit is less than the actual time 1 value, then shift the limit to reconcile:
        if (lambdaCread[m]==0 & cafeC_baset1[m]<=cafeCtarget[1]) cafelimC[m] = cafeCtarget[t] - (cafeCtarget[1] - cafeC_baset1[m]) - 0.01;
        if (lambdaTread[m]==0 & cafeT_baset1[m]<=cafeTtarget[1]) cafelimT[m] = cafeTtarget[t] - (cafeTtarget[1] - cafeT_baset1[m]) - 0.01;

        //now guess cafestat -- default is that it will stay the same as last time period (or 0 in the case of t==1)

        if (t==1 && lambdaCread[m]>0)  cafeCstat[m]=1.0; //binding according to lamdaread
        if (t==1 && lambdaTread[m]>0)  cafeTstat[m]=1.0;
        
        //other guesses (sometimes) needed to help the solver -- e.g.:
        if (t==3 && m==1) cafeCstat[m] = 0.0;
        if (t==4 && m==7) cafeTstat[m] = 1.0;
        if (t==6 && m==6) cafeTstat[m] = 1.0;
        if (t==7 && m==1) cafeCstat[m] = 1.0;
        if (t==7 && m==6) cafeCstat[m] = 1.0;

		//region-level target
        if (runmode != 1) {
            cafelimP[m] = 1.0;
            if (t==1) cafePstat[m] = 0.0;
        }
        else {
            cafelimP[m]=pavtarget[t];
            if (t==1) cafePstat[m] = 0.0;
            if (t==1 && cafeP_baset1[m]<pavtarget[t]) cafePstat[m]=1.0;
        }

	} //end loop over makes

}

void check_binding_violations(void) {
//procedure checks if there are no negative lambdas and no violated constraints.
//if it finds a negative lambda, the corresponding bindingness switch is set to 0.
//if it finds a violated constraint, the corresponding bindingness switch is set to 1.
	int m, counter;
	counter = 0;

	for (m=1;m<=NMAKE;m++) {
		if (lambdaC[m] < 0.0) {
			cafeCstat[m] = 0.0;
			lambdaC[m] = 0.0;
			counter += 1;
            if(verbose) printf("\n\n NEGATIVE LAMBDA C make %d",m);
		}
		if (lambdaT[m] < 0.0) {
			cafeTstat[m] = 0.0;
			lambdaT[m] = 0.0;
			counter += 1;
            if(verbose) printf("\n\n NEGATIVE LAMBDA T make %d",m);
		}
		if (lambdaP[m] < 0.0) {
			cafePstat[m] = 0.0;
			lambdaP[m] = 0.0;
			counter += 1;
            if(verbose) printf("\n\n NEGATIVE LAMBDA P make %d",m);
		}
	}

	for (m=1;m<=NMAKE;m++) {
		if (cafeC[m]<cafelimC[m]-0.01) {
			cafeCstat[m] = 1.0;
			counter += 1;			
            if(verbose) printf("\n\n CAFE C VIOLATED make %d",m);
		}
		if (cafeT[m]<cafelimT[m]-0.01) {
			cafeTstat[m] = 1.0;
			counter += 1;
            if(verbose) printf("\n\n CAFE T VIOLATED make %d",m);
		}
        if (cafeP[m]<cafelimP[m]-0.01) {
            cafePstat[m] = 1.0;
            counter += 1;
            if(verbose) printf("\n\n CAFE P VIOLATED make %d",m);
        }
	}
	
	if (counter > 0) binding_violations = 1;
	else binding_violations = 0;

}

void readinput(void) {
	//read input files
	FILE *infile, *lambdas, *cafetargets, *pollutionfile, *pollutionstepfile, *pollutionagefactorfile, *damages, *vmtfile, *policyfile, *tailpipecostfile, *tailpipeparamsfile, *tailpipepolicyfile, *ddafile;
	char rline[512]; //temporary string buffer for line at a time reading (lines under 512 chars)
	int m, j, g, readt, input_policy_maxt;
	double tmpbce1, tmpbce2, tmpbce3, tmpbce4, etatemp, prodemissmult;
	
	//*************** read main control file (toggles, constants, policy run controls, etc.
	infile = fopen("controlfile.txt", "r"); //open main control file
	
	fgets(rline,512,infile);
	sscanf(rline, "%s", runname); trim(runname); //string identifier for run
	fgets(rline,512,infile);
	sscanf(rline, "%s", datadir); trim(datadir); //directory to read all data from
	fgets(rline,512,infile);
	sscanf(rline, "%s", outdir); trim(outdir); //directory to write to
	fgets(rline,512,infile);
	sscanf(rline, "%s", cardatafile); trim(cardatafile); //file for car data
	fgets(rline,512,infile);
	sscanf(rline, "%s", lambdafile); trim(lambdafile); //file for lambda values
	fgets(rline,512,infile);
	sscanf(rline, "%s", elasfile); trim(elasfile); //file for elasticity data
	fgets(rline,512,infile);
	sscanf(rline, "%s", cafetargetfile); trim(cafetargetfile); //file for CAFE target data
    fgets(rline, 512, infile);
    sscanf(rline, "%s", pollutionfilename); trim(pollutionfilename); //file for pollution rates
    fgets(rline, 512, infile);
    sscanf(rline, "%s", pollutionstepfilename); trim(pollutionstepfilename); //file for pollution step multiplier (exogenous improvement for age bin 0, relative to previous time step)
    fgets(rline, 512, infile);
    sscanf(rline, "%s", pollutionagefactorfilename); trim(pollutionagefactorfilename); //file for pollution decay (relative to emissions when in age bin 0)
    fgets(rline, 512, infile);
    sscanf(rline, "%s", vmtfilename); trim(vmtfilename); //file for vmt schedule
    fgets(rline, 512, infile);
    sscanf(rline, "%s", damagesfile); trim(damagesfile); //file for damages data

	fgets(rline,512,infile);
	sscanf(rline, "%d", &runmoderead);	//0 = base, 1 = policy
    fgets(rline,512,infile);
    sscanf(rline, "%d", &nt); //number of steps to simulate
	fgets(rline,512,infile);
	sscanf(rline, "%lf", &drate); //discount rate
	fgets(rline,512,infile);
	sscanf(rline, "%ld", &ranseed);	//random seed 
	
	fgets(rline,512,infile);
	sscanf(rline, "%lf", &tmpbce1);			
	fgets(rline,512,infile);
	sscanf(rline, "%lf", &tmpbce2);			
	fgets(rline,512,infile);
	sscanf(rline, "%lf", &tmpbce3);			
	fgets(rline,512,infile);
	sscanf(rline, "%lf", &tmpbce4);			
	
    fgets(rline,512,infile);
    sscanf(rline, "%lf", &proptaxrate);
    fgets(rline,512,infile);
    sscanf(rline, "%d", &policytype);
    fgets(rline,512,infile);
    sscanf(rline, "%d", &policystart);
    fgets(rline,512,infile);
    sscanf(rline, "%d", &policyagestart);
    fgets(rline,512,infile);
    sscanf(rline, "%lf", &pollutionpolicy);
    fgets(rline,512,infile);
    sscanf(rline, "%lf", &agepolicy);
    fgets(rline, 512, infile);
    sscanf(rline, "%lf", &newcarpolicy);
    fgets(rline, 512, infile);
    sscanf(rline, "%lf", &proptaxpolicy);
    fgets(rline, 512, infile);
    sscanf(rline, "%s", tailpipecostfilename); trim(tailpipecostfilename);
    fgets(rline, 512, infile);
    sscanf(rline, "%s", tailpipeparamsfilename); trim(tailpipeparamsfilename);
    fgets(rline, 512, infile);
    sscanf(rline, "%s", tailpipepolicyfilename); trim(tailpipepolicyfilename);
    fgets(rline,512,infile);
    sscanf(rline, "%lf", &prodemissmult);
    fgets(rline,512,infile);
    sscanf(rline, "%d", &lcpolicy);
    fgets(rline, 512, infile);
    sscanf(rline, "%s", policyfilename); trim(policyfilename); //file for pollution tax policy
	fgets(rline,512,infile);
	sscanf(rline, "%d", &splitmpg);
	fgets(rline,512,infile);
    sscanf(rline, "%d", &splitused);
    fgets(rline,512,infile);
    sscanf(rline, "%d", &imperfectcomp);
    fgets(rline,512,infile);
    sscanf(rline, "%d", &dqdpusedeqm);
    fgets(rline,512,infile);
    sscanf(rline, "%lf", &tauexpectations);
    fgets(rline,512,infile);
    sscanf(rline, "%lf", &mpgexpectations);
    fgets(rline,512,infile);
    sscanf(rline, "%lf", &M_growth);
	fgets(rline,512,infile);
	sscanf(rline, "%lf", &exogmpgtech);
	fgets(rline,512,infile);
	sscanf(rline, "%lf", &gas_price); //gasoline price (expected units: $)
	fgets(rline,512,infile);
	sscanf(rline, "%lf", &theta_1);
	fgets(rline,512,infile);
	sscanf(rline, "%lf", &M_base_1); //GDP in adopting region (expected units: thousands of $)
	fgets(rline,512,infile);
	sscanf(rline, "%lf", &M_base_2); //GDP in non-adopting region (expected units: thousands of $)
	fgets(rline,512,infile);
	sscanf(rline, "%lf", &etatemp);
    fgets(rline,512,infile);
    sscanf(rline, "%d", &depreciation_mode);

	fclose(infile);

	M_1 = M_base_1;
	M_2 = M_base_2;
	M_total_1 = M_base_1;
	M_total_2 = M_base_2;

	//fill in the whole eta (scrap elasticity) vector

    for (j=1;j<=nj;j++) {
        if (age[j] <= 4 && type[j]==1 && size[j]==1) eta[j] = -0.758;
        if (age[j] >= 5 && type[j]==1 && size[j]==1) eta[j] = -0.514;
        if (age[j] <= 4 && type[j]==1 && size[j]==2) eta[j] = -0.979;
        if (age[j] >= 5 && type[j]==1 && size[j]==2) eta[j] = -0.500;
        if (age[j] <= 4 && type[j]==2 && size[j]==1) eta[j] = -0.816;
        if (age[j] >= 5 && type[j]==2 && size[j]==1) eta[j] = -0.811;
        if (age[j] <= 4 && type[j]==2 && size[j]==2) eta[j] = -0.617;
        if (age[j] >= 5 && type[j]==2 && size[j]==2) eta[j] = -1.018;
    }
    
    for (j=1;j<=nj;j++) eta[j] = etatemp * eta[j];
    
    //fill in the production damages vector in thousands of dollars (600 from IO table approach, expressed in $2019)
    for (j=1;j<=nj;j++) {
        if (j<=njnew) proddamage[j] = 600.0 * prodemissmult;
        else proddamage[j] = 0.0;
    }

	//expand the bce parameters into all new cars
	for (j=1;j<=njnew;j++) {
		if (type[j]==1 && size[j]==1) bce[j] = tmpbce1;
		if (type[j]==1 && size[j]==2) bce[j] = tmpbce2;
		if (type[j]==2 && size[j]==1) bce[j] = tmpbce3;
        if (type[j]==2 && size[j]==2) bce[j] = tmpbce4;
	}
    
    //record the base run name (if base mode)
    if(runmoderead==0) sprintf(baserunname,"%s",runname);
	
	
	// read file with the starting lambdas for the producer problem in it
	sprintf(sbuf,"%s%s",datadir,lambdafile);
	lambdas = fopen(sbuf, "r");
    
	for (m=1; m<=NMAKE; m++) {
		fscanf(lambdas, "%lf %lf", &lambdaCread[m], &lambdaTread[m]);
	}
	fclose(lambdas);	
	
	for (j=1; j<=nt; j++) {
        pavtarget[j] = 1.0; //unused, can read this to set alternative targets over regions
	}

	sprintf(sbuf,"%s%s",datadir,cafetargetfile);
	cafetargets = fopen(sbuf, "r"); 
	
	for (j=1; j<=nt; j++) {
		fscanf(cafetargets, "%lf %lf", &cafeCtarget[j], &cafeTtarget[j]);
	}
	fclose(cafetargets);	

	sprintf(sbuf, "%s%s", datadir, damagesfile);
	damages = fopen(sbuf, "r");

	fgets(rline, 512, damages);
	sscanf(rline, "%lf", &damage_co); // CO damages in $/metric ton
	fgets(rline, 512, damages);
	sscanf(rline, "%lf", &damage_hc); // HC damages in $/metric ton
	fgets(rline, 512, damages);
	sscanf(rline, "%lf", &damage_nox); // NOx damages in $/metric ton
	fclose(damages);

    //read in pollution data for the base period -- pollution from each vehicle per mile in the base time step
    sprintf(sbuf, "%s%s", datadir, pollutionfilename);
    pollutionfile = fopen(sbuf, "r");
    
    fgets(rline, 512, pollutionfile);
    sscanf(rline, "%d", &input_phi_maxt); // number of steps to expect in the file
    
    for (j = 1; j <= nj; j++) {
        for (readt=1; readt<=input_phi_maxt; readt++) {
            fscanf(pollutionfile, "%lf %lf %lf", &jt_phi_co[j][readt], &jt_phi_hc[j][readt], &jt_phi_nox[j][readt]);
        }
    }
    fclose(pollutionfile);
    
    //read in pollution step data (exogenous, new vehicles only)
    sprintf(sbuf, "%s%s", datadir, pollutionstepfilename);
    pollutionstepfile = fopen(sbuf, "r");
    
    for (j = 1; j <= njnew; j++) {
        fscanf(pollutionstepfile, "%lf %lf %lf", &phi_co_newcarstep[j], &phi_hc_newcarstep[j], &phi_nox_newcarstep[j]);
    }
    fclose(pollutionstepfile);
    
    
    //read in pollution age factor data (emissions relative to emissions when vehicle was in age bin 0)
    sprintf(sbuf, "%s%s", datadir, pollutionagefactorfilename);
    pollutionagefactorfile = fopen(sbuf, "r");
    
    for (j = 1; j <= nj; j++) {
        fscanf(pollutionagefactorfile, "%lf %lf %lf", &phi_co_agefactor[j], &phi_hc_agefactor[j], &phi_nox_agefactor[j]);
    }
    fclose(pollutionagefactorfile);
    

	//read in VMT by age-model schedule
    sprintf(sbuf, "%s%s", datadir, vmtfilename);
    vmtfile = fopen(sbuf, "r");
    
	for (j = 1; j <= nj; j++) {
		fscanf(vmtfile, "%lf", &vmt[j]);
	}
	fclose(vmtfile);
    
    //read in pollution policy (tau) direct;y
    if (policytype == 10) {
        sprintf(sbuf, "%s%s", datadir, policyfilename);
        policyfile = fopen(sbuf, "r");
        
        input_policy_maxt = 50;  //50 time steps in file as default (values past nt will be ignored)
        
        for (j = 1; j <= nj; j++) {
            for (readt=1; readt<=input_policy_maxt; readt++) {
                fscanf(policyfile, "%lf", &tau_read[j][readt]);
            }
        }
        fclose(policyfile);
    }
    
    //read in baseline tailpipe cost file (new cars)
    sprintf(sbuf, "%s%s", datadir, tailpipecostfilename);
    tailpipecostfile = fopen(sbuf, "r");
    
    input_policy_maxt = 30;
    for (j = 1; j <= njnew; j++) {
        for (readt=1; readt<=input_policy_maxt; readt++) {
            fscanf(tailpipecostfile, "%lf", &basetailpipecost[j][readt]);
        }
    }
    
    //extend to max t considered (only relevant if nt is set very high)
    for (readt=input_policy_maxt+1; readt<=nt+NAGE+1; readt++) {
        basetailpipecost[j][readt] = basetailpipecost[j][input_policy_maxt];
    }
    
    fclose(tailpipecostfile);
    
    
    //read in tailpipe calibrated cost params
    sprintf(sbuf, "%s%s", datadir, tailpipeparamsfilename);
    tailpipeparamsfile = fopen(sbuf, "r");
    
    fscanf(tailpipeparamsfile, "%lf", &tailpipez);
    
    for (j = 1; j <= njnew; j++) {
        fscanf(tailpipeparamsfile, "%lf", &tailpipev[j]);
    }
    
    fclose(tailpipeparamsfile);
    
    
    //read in tailpipe policy
    sprintf(sbuf, "%s%s", datadir, tailpipepolicyfilename);
    tailpipepolicyfile = fopen(sbuf, "r");
    
    input_policy_maxt = 50;  //50 time steps in file as default (will be held constant past that as needed)
    
    for (j = 1; j <= njnew; j++) {
        for (readt=1; readt<=input_policy_maxt; readt++) {
            fscanf(tailpipepolicyfile, "%lf", &tailpipepolicy[j][readt]);
        }
    }
    //extend to max t considered (only relevant if nt is set very high)
    for (readt=input_policy_maxt+1; readt<=nt+NAGE+1; readt++) {
        tailpipepolicy[j][readt] = tailpipepolicy[j][input_policy_maxt];
    }
    
    fclose(tailpipepolicyfile);
    
    sprintf(sbuf, "%s%s", datadir, "ddafile.txt");
    ddafile = fopen(sbuf, "r");

    
    for (j = 0; j <= nj; j++) {
        for (g=1; g<=ng; g++) {
            fscanf(ddafile, "%lf", &dda_base[j][g]);
        }
    }
    fclose(ddafile);
    
    
}

void definecar(void) {
    // read car data file, initialize any other car indexed variables
    // calibrate the b parameters in the scrap function
    // calibrate the pollution model
    // impute rj vector for the baseline
    // impute payments for profit model in the baseline
    
    int j, xt, jx, ax, jxx, jnew;
    FILE *carfile;
    double dum;
    double rj0[nj+1];
    double pj0[nj+1];
    double hj0[nj+1];
    double pscrap0[nj+1];
    double mpg0[nj+1];
    double milesremaining0[nj+1];
    double dvalde0[nj+1];
    double xpj0[nj+1], xhj0[nj+1], xpscrap0[nj+1], xmilesremaining0[nj+1];
    double temppaytotal1[nj+1], temppaytotal2[nj+1];
    double jt_weightedphi[njnew+1][nt+NAGE+2], jt_weightedphibase[njnew+1][nt+NAGE+2];
    
    
    sprintf(sbuf,"%s%s",datadir,cardatafile);
    carfile = fopen(sbuf, "r"); //open the car chars input file
    
    for (j=1; j<=nj; j++) {
        fscanf(carfile, "%lf %lf %lf %lf %lf %lf", &dum, &mpg0[j], &pj0[j], &pscrap0[j], &based1[j], &based2[j]);
        pj0[j] = pj0[j] / 1000.0;
        
        if (j>njnew && pj0[j]-pj0[j-njnew] > 0) {
            printf("\n\nInconsistent input price data, price increase at %d",j);
            pj0[j]=pj0[j-njnew];
        }
        
    }
    
    fclose(carfile);
    
    //calibrate the b[j] constant terms in the scrap function to baseline prices and scrap rates:
    //(solves the scrap function for b)
    for (j=njnew+1;j<=nj;j++) {
        b[j] = pscrap0[j] / pow(pj0[j],eta[j]);
    }
        
    // calculate all the pollution "phi" values for all time steps (including projecting out past the end of the run, because we need those for implicit rental price calculations when tauexpectations>0)
    
    // assign implicit emissions when new for old cars in the base year
    for (j=1; j<=nj; j++) {
        jt_phi_co_when_new[j][1] = jt_phi_co[j][1]/phi_co_agefactor[j];
        jt_phi_hc_when_new[j][1] = jt_phi_hc[j][1]/phi_hc_agefactor[j];
        jt_phi_nox_when_new[j][1] = jt_phi_nox[j][1]/phi_nox_agefactor[j];
    }
    
    for (xt=2; xt<=nt+NAGE+1; xt++) {
        //assign emissions for vehicles in future fleets
        
        //assign when_new emissions to successive ages through time, keeping them with the vintage
        for (j = njnew+1; j <= nj; j++) {
            jt_phi_co_when_new[j][xt] = jt_phi_co_when_new[j-njnew][xt-1];
            jt_phi_hc_when_new[j][xt] = jt_phi_hc_when_new[j-njnew][xt-1];
            jt_phi_nox_when_new[j][xt] = jt_phi_nox_when_new[j-njnew][xt-1];
        }
        if(xt<=input_phi_maxt) {
            //use emissions as read for future new vehicles
            for (j = 1; j <= njnew; j++) {
                jt_phi_co_when_new[j][xt] = jt_phi_co[j][xt];
                jt_phi_hc_when_new[j][xt] = jt_phi_hc[j][xt];
                jt_phi_nox_when_new[j][xt] = jt_phi_nox[j][xt];
            }
            
        }
        else {
            //compute emissions for future new vehicles that are beyond input_phi_maxt
            for (j = 1; j <= njnew; j++) {
                jt_phi_co_when_new[j][xt] = jt_phi_co_when_new[j][xt-1] * phi_co_newcarstep[j];
                jt_phi_hc_when_new[j][xt] = jt_phi_hc_when_new[j][xt-1] * phi_hc_newcarstep[j];
                jt_phi_nox_when_new[j][xt] = jt_phi_nox_when_new[j][xt-1] * phi_nox_newcarstep[j];
            }
            //compute all emissions (note agefactor = 1 for new vehicles)
            for (j = 1; j <= nj; j++) {
                jt_phi_co[j][xt] = jt_phi_co_when_new[j][xt] * phi_co_agefactor[j];
                jt_phi_hc[j][xt] = jt_phi_hc_when_new[j][xt] * phi_hc_agefactor[j];
                jt_phi_nox[j][xt] = jt_phi_nox_when_new[j][xt] * phi_nox_agefactor[j];
            }
        }
    }
    
    // calibrate the cost functions for tailpipe emissions improvements (phis are all at baseline here)
    
    for (j = 1; j <= njnew; j++) {
        for (xt=1; xt<=nt+NAGE+1; xt++) {
            jt_weightedphibase[j][xt] = damage_co * jt_phi_co[j][xt] + damage_hc * jt_phi_hc[j][xt] + damage_nox * jt_phi_nox[j][xt];
            tailpipew[j][xt] = basetailpipecost[j][xt] - (   pow(tailpipez,xt) * tailpipev[j] * (jt_weightedphibase[j][1]/jt_weightedphibase[j][xt] - 1)   );
        }
    }

    // apply tailpipe policy (if any) to re-compute emissions for all vehicles
    for (jnew = 1; jnew <= njnew; jnew++) {
        for (xt=1; xt<=nt+NAGE+1; xt++) {
            // follow this vintage as it ages
            for (ax=0; ax<=18; ax++) {
                if (ax > nt+NAGE+1-xt) continue;
                jt_phi_co[jnew+ax*njnew][xt+ax] *= 1 - tailpipepolicy[jnew][xt];
                jt_phi_hc[jnew+ax*njnew][xt+ax] *= 1 - tailpipepolicy[jnew][xt];
                jt_phi_nox[jnew+ax*njnew][xt+ax] *= 1 - tailpipepolicy[jnew][xt];
                jt_phi_co_when_new[jnew+ax*njnew][xt+ax] *= 1 - tailpipepolicy[jnew][xt];
                jt_phi_hc_when_new[jnew+ax*njnew][xt+ax] *= 1 - tailpipepolicy[jnew][xt];
                jt_phi_nox_when_new[jnew+ax*njnew][xt+ax] *= 1 - tailpipepolicy[jnew][xt];
            }
        }
    }

    
    // compute tailpipe cost (reproduces baseline if tailpipepolicy = 0)
    for (j = 1; j <= njnew; j++) {
        for (xt=1; xt<=nt+NAGE+1; xt++) {
            jt_weightedphi[j][xt] = damage_co * jt_phi_co[j][xt] + damage_hc * jt_phi_hc[j][xt] + damage_nox * jt_phi_nox[j][xt];
            jt_tailpipecost[j][xt] = pow(tailpipez,xt) * tailpipev[j] * (jt_weightedphibase[j][1]/jt_weightedphi[j][xt] - 1) + tailpipew[j][xt];
        }
    }

    //initialize car pollution data and tailpipe costs for time period 1 (pulling from the jt matrix)
    for (j=1; j<=nj; j++) {
        phi_co[j] = jt_phi_co[j][1];
        phi_hc[j] = jt_phi_hc[j][1];
        phi_nox[j] = jt_phi_nox[j][1];
        phi_co_when_new[j] = jt_phi_co_when_new[j][1];
        phi_hc_when_new[j] = jt_phi_hc_when_new[j][1];
        phi_nox_when_new[j] = jt_phi_nox_when_new[j][1];
    }
    for (j = 1; j <= njnew; j++) {
        tailpipecost[j] = jt_tailpipecost[j][1];
    }
    
    
    // calibrate vehicle prices at t = 0 (baseline tailpipe cost is relative to t=1, so 0)
    for (j=nj;j>=1;j--) {
        
        //loop through from future oldest "jx" versions of car j backward
        for (ax=18; ax>=age[j]; ax--) {
            jx = j+(ax-age[j])*njnew;
            if (ax==18) {
                if (j==jx) {
                    exprj1[jx][1] = pj0[jx];
                }
                else {
                    exprj1[jx][1+ax-age[j]] = exprj1[jx][1] + mpgexpectations*0.001*vmt[jx]*gas_price * (1/mpg0[jx] - 1/mpg0[j]);
                }
                xpj0[jx] = exprj1[jx][1+ax-age[j]]; //this is the price we expect car j to have when it is 18 years old (and so has become car "jx")
                xmilesremaining0[jx] = vmt[jx];
            }
            if (ax>0 && ax<18) {
                if (j==jx) {
                    exprj1[jx][1] = pj0[jx] - (1.0 - xpscrap0[jx+njnew])*(xpj0[jx+njnew]-xhj0[jx+njnew])/(1.0 + drate);
                }
                else {
                    exprj1[jx][1+ax-age[j]] = exprj1[jx][1] + mpgexpectations*0.001*vmt[jx]*gas_price * (1/mpg0[jx] - 1/mpg0[j]);
                }
                xpj0[jx] = exprj1[jx][1+ax-age[j]] + (1.0 - xpscrap0[jx+njnew])*(xpj0[jx+njnew]-xhj0[jx+njnew])/(1.0 + drate);
                xmilesremaining0[jx] = vmt[jx] + (1.0 - xpscrap0[jx+njnew]) * xmilesremaining0[jx+njnew] / (1.0 + drate);
                
            }
            
            //new cars
            if (ax==0) {
                //determine rj for new cars -- this point in the loop is only reached for currently-new cars; jx=j here
                
                rj0[jx] = pj0[jx] - (1.0 - xpscrap0[jx+njnew]) * (xpj0[jx+njnew]-xhj0[jx+njnew])/(1.0 + drate);
                
                milesremaining0[jx] = vmt[jx] + (1.0 - xpscrap0[jx+njnew]) * xmilesremaining0[jx+njnew] / (1.0 + drate);
                
                pscrap0[jx] = 0.0;
                hj0[jx] = 0.0;
                
                fracremain0[j] = 1.0;  //fraction of each new car expected to remain over time (age in this notation)
                for (jxx=j+njnew; jxx<=nj; jxx+=njnew) {
                    fracremain0[jxx] = fracremain0[jxx-njnew] * (1-xpscrap0[jxx]);
                }
            }
            
            //determine scrap rate
            if (ax>0) {
                if (xpj0[jx]>0) {
                    xpscrap0[jx] = b[jx] * pow(xpj0[jx],eta[jx]);
                    xhj0[jx] = (-1.0 * pow(b[jx],-1/eta[jx]) * eta[jx] + b[jx] * eta[jx] * pow(xpj0[jx],1+eta[jx])) /
                    ((1+eta[jx]) * (-1.0 + b[jx] * pow(xpj0[jx],eta[jx])));
                }
                else {
                    xpscrap0[jx] = 1.0;
                    xhj0[jx] = 0.0;
                }
                xpscrap0[jx] = DMIN(xpscrap0[jx],0.9999);
            }
            
            //if jx=j then we've finished working our way backward through future ages and can assign values to the current-period arrays
            if (jx==j && ax>0) {
                rj0[j] = exprj1[j][1];
                pscrap0[j] = xpscrap0[j];
                hj0[j] = xhj0[j];
                milesremaining0[j] = xmilesremaining0[j];
            }
            if (jx==j && ax<18) {
                //save current expectations about prices etc. one year forward (for use in capital gains calculations)
                //also assumed to be same from last (i.e. t=0) period in steady state
                lastexppj1[j+njnew] = exppj1[j+njnew] = xpj0[j+njnew];
                lastexppscrap1[j+njnew] = exppscrap1[j+njnew] = xpscrap0[j+njnew];
                lastexphj1[j+njnew] =  exphj1[j+njnew] = xhj0[j+njnew];
                lastexppj2[j+njnew] = exppj2[j+njnew] = xpj0[j+njnew];
                lastexppscrap2[j+njnew] = exppscrap2[j+njnew] = xpscrap0[j+njnew];
                lastexphj2[j+njnew] = exphj2[j+njnew] = xhj0[j+njnew];
            }
        }
    }
    
    // dvalde0 for new cars at mpg0
    for (j=1; j<=njnew; j++) dvalde0[j] = milesremaining0[j] * 0.001 * pow(mpg0[j],-2) * gas_price;
    
    //populate past payments
    //(these were determined when existing stock was purchased new, so will be deducted from income regardless of policy choice)
    for (j=1; j<=nj; j++) {
        //payment schedule at baseline scrappage and rental prices
        payment1[j][0] = fracremain0[j] * ( rj0[j] - hj0[j] +
                                           mpgexpectations*0.001*vmt[j]*gas_price*(1/mpg0[j]-1/mpg0[j-age[j]*njnew]) );
    }
    // reconcile the (typically very small) differences coming from endogeneity of scrap rates and mpg adjustments
    for (j=1; j<=nj; j++) {
        if (age[j]==0) temppaytotal1[j] = payment1[j][0];
        else temppaytotal1[j-njnew*age[j]] += payment1[j][0] / pow(1.0+drate,age[j]);
    }
    for (j=1; j<=nj; j++) {
        payment1[j][0] = payment1[j][0] * pj0[j-njnew*age[j]]/temppaytotal1[j-njnew*age[j]] ;
    }
    //now populate the array of payments on cars purchased before t=1 (payments for cars purchased from t=1 onward will populate as the model runs)
    for (j=1; j<=nj; j++) {
        for (xt=1; xt<=age[j]; xt++) {
            payment1[j][xt] = payment1[j][0];
            payment2[j][xt] = payment1[j][xt];
        }
    }

    //store implied demand when new for all baseline cars
    for (j=1; j<=nj; j++) {
        d1_when_new[j] = based1[j] / fracremain0[j];
        d2_when_new[j] = based2[j] / fracremain0[j];
        lastd1[j]=based1[j];
        lastd2[j]=based2[j];
    }
    
    
    //set two regions equal to start with
    for (j=1; j<=nj; j++) {
        baserj1[j] = rj0[j];
        baserj2[j] = rj0[j];
        basepj1[j] = pj0[j];
        basepj2[j] = pj0[j];
        basepscrap1[j] = pscrap0[j];
        basepscrap2[j] = pscrap0[j];
        basehj1[j] = hj0[j];
        basehj2[j] = hj0[j];
        mpg1[j] = mpg0[j];
        mpg2[j] = mpg0[j];
        dvalde1[j] = dvalde0[j];
        dvalde2[j] = dvalde0[j];
        milesremaining1[j] = milesremaining0[j];
        milesremaining2[j] = milesremaining0[j];
        hj1[j] = basehj1[j];
        hj2[j] = basehj2[j];
        rj1[j] = baserj1[j];
        rj2[j] = baserj2[j];
        pj1[j] = basepj1[j];
        pj2[j] = basepj2[j];
        pscrap1[j] = basepscrap1[j];
        pscrap2[j] = basepscrap2[j];
    }
    
}


void definecalib(void) {
    //calibration program to assign rho and alphas (using baseline prices and quantities)
    //also reads policy utilities (if relevant, for use in EV runs)
    
	int dim1;
	int dim2;
	int dim3;
	int dim4;
	int j;
    double q_temp_1;
    double q_temp_2;
    double mpg4dim[NAGE + 1][NTYPE + 1][NSIZE + 1][NMAKE + 1];
    double vmt4dim[NAGE + 1][NTYPE + 1][NSIZE + 1][NMAKE + 1];
    double proptax4dim[NAGE + 1][NTYPE + 1][NSIZE + 1][NMAKE + 1];
    double r_0_1[NAGE + 1][NTYPE + 1][NSIZE + 1][NMAKE + 1];
    double r_0_2[NAGE + 1][NTYPE + 1][NSIZE + 1][NMAKE + 1];
    double q_0_1[NAGE + 1][NTYPE + 1][NSIZE + 1][NMAKE + 1];
    double q_0_2[NAGE + 1][NTYPE + 1][NSIZE + 1][NMAKE + 1];
    double tap_0_1[NAGE + 1][NTYPE + 1][NSIZE + 1][NMAKE + 1];
    double tap_0_2[NAGE + 1][NTYPE + 1][NSIZE + 1][NMAKE + 1];
    double q_comp_1_1[NAGE + 1][NTYPE + 1][NSIZE + 1];
    double q_comp_1_2[NAGE + 1][NTYPE + 1][NSIZE + 1];
    double q_comp_2_1[NTYPE + 1][NSIZE + 1];
    double q_comp_2_2[NTYPE + 1][NSIZE + 1];
    double q_comp_3_1[NTYPE + 1];
    double q_comp_3_2[NTYPE + 1];
    double q_comp_4_1;
    double q_comp_4_2;
    double q_comp_5_1;
    double q_comp_5_2;
    int g;
    double totddaj[nj+1];
    double totdda, totmodel;
    double mvecguess[nj+ng+1], mvecdiffs[nj+ng+1];
	FILE *elasticityfile;
	FILE *bfile;
	
	sprintf(sbuf,"%s%s",datadir,elasfile);
	elasticityfile = fopen(sbuf, "r"); //open the elasticity input file

    // Read the elasticity parameters for nest 1 (for all t,s,a)
    for (dim1=1; dim1<=19; dim1++) {
        for (dim2=1; dim2<=2; dim2++) {
            for (dim3=1; dim3<=2; dim3++) {
                fscanf(elasticityfile, "%lf", &rho_1[dim1][dim2][dim3]);
            }
        }
    }
    
    // Read the elasticity parameters for nest 2 (for all t,s)
    for (dim2=1; dim2<=2; dim2++) {
        for (dim3=1; dim3<=2; dim3++) {
            fscanf(elasticityfile, "%lf", &rho_2[dim2][dim3]);
        }
    }
    
    // Read the elasticity parameters for nest 3 (for all t)
    for (dim2=1; dim2<=2; dim2++) {
        fscanf(elasticityfile, "%lf", &rho_3[dim2]);
    }
    
    // Read the elasticity parameter for nest 4 (rho_v in demand writeup)
    fscanf(elasticityfile, "%lf", &rho_4);
    // Read the elasticity parameter for nest 5 (rho_u in demand writeup))
    fscanf(elasticityfile, "%lf", &rho_5);
    
    fclose(elasticityfile);

    
    
    // put the VMT, MPG, rj, and dj data into the 4-dimensional array structure needed for the calibration program:
    j=1;
    for (dim1 = 1; dim1 <= 19; dim1++) {
        for (dim2 = 1; dim2 <= 2; dim2++) {
            for (dim3 = 1; dim3 <= 2; dim3++) {
                for (dim4 = 1; dim4 <= 7; dim4++) {
                    vmt4dim[dim1][dim2][dim3][dim4] = vmt[j];
                    mpg4dim[dim1][dim2][dim3][dim4] = mpg1[j]; // uniform MPG over regions for baseline calib
                    r_0_1[dim1][dim2][dim3][dim4] = baserj1[j];
                    r_0_2[dim1][dim2][dim3][dim4] = baserj2[j];
                    q_0_1[dim1][dim2][dim3][dim4] = based1[j];
                    q_0_2[dim1][dim2][dim3][dim4] = based2[j];
                    proptax4dim[dim1][dim2][dim3][dim4] = 2 * proptaxrate * basepj1[j] * 1000;
                    j++;
                }
            }
        }
    }
    

    // Compute composite vehicle quantity for lowest nest (nest 1)
    for (dim2=1; dim2<=2; dim2++) {
        for (dim3=1; dim3<=2; dim3++) {
            for (dim1=1; dim1<=19; dim1++) {
                q_temp_1 = 0;
                q_temp_2 = 0;
                for (dim4=1; dim4<=7; dim4++) {
                    tap_0_1[dim1][dim2][dim3][dim4] = r_0_1[dim1][dim2][dim3][dim4]
                    + (vmt4dim[dim1][dim2][dim3][dim4] / mpg4dim[dim1][dim2][dim3][dim4]) * gas_price / 1000.0
                    + proptax4dim[dim1][dim2][dim3][dim4]/ 1000.0;
                    tap_0_2[dim1][dim2][dim3][dim4] = r_0_2[dim1][dim2][dim3][dim4]
                    + (vmt4dim[dim1][dim2][dim3][dim4] / mpg4dim[dim1][dim2][dim3][dim4]) * gas_price / 1000.0
                    + proptax4dim[dim1][dim2][dim3][dim4]/ 1000.0;
                    q_temp_1 += q_0_1[dim1][dim2][dim3][dim4] * tap_0_1[dim1][dim2][dim3][dim4] / tap_comp_1_1[dim1][dim2][dim3];
                    q_temp_2 += q_0_2[dim1][dim2][dim3][dim4] * tap_0_2[dim1][dim2][dim3][dim4] / tap_comp_1_2[dim1][dim2][dim3];
                }
                q_comp_1_1[dim1][dim2][dim3] = q_temp_1;
                q_comp_1_2[dim1][dim2][dim3] = q_temp_2;
            }
        }
    }
    
    // Compute distribution parameters for lowest nest (nest 1)
    for (dim2=1; dim2<=2; dim2++) {
        for (dim3=1; dim3<=2; dim3++) {
            for (dim1=1; dim1<=19; dim1++) {
                for (dim4=1; dim4<=7; dim4++) {
                    alpha_1_1[dim1][dim2][dim3][dim4] = (tap_0_1[dim1][dim2][dim3][dim4] / tap_comp_1_1[dim1][dim2][dim3]) * pow((q_0_1[dim1][dim2][dim3][dim4] / q_comp_1_1[dim1][dim2][dim3]),(1 - rho_1[dim1][dim2][dim3]));
                    alpha_1_2[dim1][dim2][dim3][dim4] = (tap_0_2[dim1][dim2][dim3][dim4] / tap_comp_1_2[dim1][dim2][dim3]) * pow((q_0_2[dim1][dim2][dim3][dim4] / q_comp_1_2[dim1][dim2][dim3]),(1 - rho_1[dim1][dim2][dim3]));
                }
            }
        }
    }
    
    // Compute composite vehicle quantity for nest 2
    for (dim2=1; dim2<=2; dim2++) {
        for (dim3=1; dim3<=2; dim3++) {
            q_temp_1 = 0;
            q_temp_2 = 0;
            for (dim1=1; dim1<=19; dim1++) {
                q_temp_1 += q_comp_1_1[dim1][dim2][dim3]*tap_comp_1_1[dim1][dim2][dim3]/tap_comp_2_1[dim2][dim3];
                q_temp_2 += q_comp_1_2[dim1][dim2][dim3]*tap_comp_1_2[dim1][dim2][dim3]/tap_comp_2_2[dim2][dim3];
            }
            q_comp_2_1[dim2][dim3] = q_temp_1;
            q_comp_2_2[dim2][dim3] = q_temp_2;
        }
    }
    
    // Compute distribution parameters for nest 2
    for (dim2=1; dim2<=2; dim2++) {
        for (dim3=1; dim3<=2; dim3++) {
            for (dim1=1; dim1<=19; dim1++) {
                alpha_2_1[dim1][dim2][dim3] = (tap_comp_1_1[dim1][dim2][dim3] / tap_comp_2_1[dim2][dim3]) * pow((q_comp_1_1[dim1][dim2][dim3] / q_comp_2_1[dim2][dim3]),(1 - rho_2[dim2][dim3]));
                alpha_2_2[dim1][dim2][dim3] = (tap_comp_1_2[dim1][dim2][dim3] / tap_comp_2_2[dim2][dim3]) * pow((q_comp_1_2[dim1][dim2][dim3] / q_comp_2_2[dim2][dim3]),(1 - rho_2[dim2][dim3]));
            }
        }
    }
    
    // Compute composite vehicle quantity for nest 3
    for (dim2=1; dim2<=2; dim2++) {
        q_temp_1 = 0;
        q_temp_2 = 0;
        for (dim3=1; dim3<=2; dim3++) {
            q_temp_1 += q_comp_2_1[dim2][dim3]*tap_comp_2_1[dim2][dim3]/tap_comp_3_1[dim2];
            q_temp_2 += q_comp_2_2[dim2][dim3]*tap_comp_2_2[dim2][dim3]/tap_comp_3_2[dim2];
        }
        q_comp_3_1[dim2] = q_temp_1;
        q_comp_3_2[dim2] = q_temp_2;
    }
    
    // Compute distribution parameters for nest 3
    for (dim2=1; dim2<=2; dim2++) {
        for (dim3=1; dim3<=2; dim3++) {
            alpha_3_1[dim2][dim3] = (tap_comp_2_1[dim2][dim3] / tap_comp_3_1[dim2]) * pow((q_comp_2_1[dim2][dim3] / q_comp_3_1[dim2]),(1 - rho_3[dim2]));
            alpha_3_2[dim2][dim3] = (tap_comp_2_2[dim2][dim3] / tap_comp_3_2[dim2]) * pow((q_comp_2_2[dim2][dim3] / q_comp_3_2[dim2]),(1 - rho_3[dim2]));
        }
    }
    
    // Compute composite vehicle quantity for nest 4
    q_temp_1 = 0;
    q_temp_2 = 0;
    for (dim2=1; dim2<=2; dim2++) {
        q_temp_1 += q_comp_3_1[dim2]*tap_comp_3_1[dim2]/tap_comp_4_1;
        q_temp_2 += q_comp_3_2[dim2]*tap_comp_3_2[dim2]/tap_comp_4_2;
    }
    q_comp_4_1 = q_temp_1;
    q_comp_4_2 = q_temp_2;
    
    // Compute distribution parameters for nest 4
    for (dim2=1; dim2<=2; dim2++) {
        alpha_4_1[dim2] = (tap_comp_3_1[dim2] / tap_comp_4_1) * pow((q_comp_3_1[dim2] / q_comp_4_1),(1 - rho_4));
        alpha_4_2[dim2] = (tap_comp_3_2[dim2] / tap_comp_4_2) * pow((q_comp_3_2[dim2] / q_comp_4_2),(1 - rho_4));
    }
    
    // Compute consumption good quantity for nest 5
    q_comp_5_1 = M_1;
    q_comp_5_2 = M_2;
    
    // Compute distribution parameters for nest 5
    alpha_5_1[1] = (tap_comp_4_1 / tap_comp_5_1) * pow((q_comp_4_1 / q_comp_5_1),(1 - rho_5));
    alpha_5_1[2] = (tap_x_1 / tap_comp_5_1) * pow((((tap_comp_5_1*q_comp_5_1 - tap_comp_4_1*q_comp_4_1)/tap_x_1) / q_comp_5_1),(1 - rho_5));
    alpha_5_2[1] = (tap_comp_4_2 / tap_comp_5_2) * pow((q_comp_4_2 / q_comp_5_2),(1 - rho_5));
    alpha_5_2[2] = (tap_x_2 / tap_comp_5_2) * pow((((tap_comp_5_2*q_comp_5_2 - tap_comp_4_2*q_comp_4_2)/tap_x_2) / q_comp_5_2),(1 - rho_5));
 
    
    //rescale the baseline dda values such that they match baseline quantities (allows baseline dda file to come from another dataset)
    //first get total cars in the dda file and compare to model to scale everything together
    totmodel = 0.0;
    totdda = 0.0;
    for (j = 1; j <= nj; j++) {
        totmodel += based1[j]+based2[j];
        for (g = 1; g <= ng; g++) {
            totdda += dda_base[j][g];
        }
    }
    for (j = 0; j <= nj; j++) {
        for (g = 1; g <= ng; g++) dda_base[j][g] = dda_base[j][g] * totmodel/totdda;
    }
    
    
    // de-aggregate vehicles by income (holding shares proportional to baseline time 0)
    // note: the model here is of a representative agent, so there are many ways to de-aggregate, this represents a simple possibility that preserves vehicle totals and proportionality to the disaggregation (from data) in the baseline
    // target vectors:
    
    for (j = 1; j <= nj; j++) {
        dda_jtarg[j] = based1[j] + based2[j];
    }
    
    // group sizes
    for (g=1; g <= ng; g++) {
        dda_gtarg[g] = 0.0;
        for (j = 0; j <= nj; j++) dda_gtarg[g] += dda_base[j][g];
    }
    
    for (j=1; j<=nj+ng+1; j++) {
        mvecguess[j]=1.0;
    }
    broydn2(mvecguess,nj+ng,&retvaldda,dda_objective);
    if (retvaldda!=0) printf("\nERROR in income deaggregation at initial baseline");
    dda_objective(nj+ng,mvecguess,mvecdiffs);

    //now overwrite the base with the calibrated one
    for (j = 0; j <= nj; j++) {
        for (g = 1; g <= ng; g++) dda_base[j][g] = dda[j][g];
    }
}


void initialization(void) {
	//set composite price to 1; specify income, VMT, gas price and scrap elasticity
	int dim1;
	int dim2;
	int dim3;
	int dim4;
    int j;
	
	//fill make[j], age[j], size[j] and type[j]
	for (dim1=1; dim1<=19; dim1++) {
		for (dim2=1; dim2<=2; dim2++) {
			for (dim3=1; dim3<=2; dim3++) {
				for (dim4=1; dim4<=7; dim4++) {
					type[((dim1 - 1)*(NTYPE * NSIZE * NMAKE) + (dim2 - 1)*(NSIZE * NMAKE) + (dim3 - 1)*NMAKE + dim4)] = dim2;
					size[((dim1 - 1)*(NTYPE * NSIZE * NMAKE) + (dim2 - 1)*(NSIZE * NMAKE) + (dim3 - 1)*NMAKE + dim4)] = dim3;
					age[((dim1 - 1)*(NTYPE * NSIZE * NMAKE) + (dim2 - 1)*(NSIZE * NMAKE) + (dim3 - 1)*NMAKE + dim4)] = dim1 - 1;
					make[((dim1 - 1)*(NTYPE * NSIZE * NMAKE) + (dim2 - 1)*(NSIZE * NMAKE) + (dim3 - 1)*NMAKE + dim4)] = dim4;
				}
			}
		}
	}
	
	//fill mc[make][vehicle] (1= sm car, 2= lg car, 3= sm truck, 4= lg truck)
	for (j=1;j<=njnew;j++) {
		mc[make[j]][(type[j]-1)*2+size[j]] = j;
	}	
	
	//now all composite prices are set to 1
	for (dim2=1; dim2<=2; dim2++) {
		for (dim3=1; dim3<=2; dim3++) {
			for (dim1=1; dim1<=19; dim1++) {
				tap_comp_1_1[dim1][dim2][dim3] = 1.0;
				tap_comp_1_2[dim1][dim2][dim3] = 1.0;
				tap_comp_2_1[dim2][dim3] = 1.0;
				tap_comp_2_2[dim2][dim3] = 1.0;
				tap_comp_3_1[dim2] = 1.0;
				tap_comp_3_2[dim2] = 1.0;
			}
		}
	}
	
	tap_comp_4_1 = 1.0;
	tap_comp_4_2 = 1.0;
	tap_comp_5_1 = 1.0;
	tap_comp_5_2 = 1.0;
	
	//normalize the price of "other goods" to 1 (in both regions)
	tap_x_1 = 1.0;
	tap_x_2 = 1.0;
	

	//set commonmpg[j] to zero initially
	for (j=1;j<=nj;j++) commonmpg[j] = 0.0;

	t = 0;
	profit1 = 0;
	profit2 = 0;
	baseprofit1 = 0;
	baseprofit2 = 0;

	dqdrnew1 = dmatrix(1,nj,1,nj);
	dqdrnew2 = dmatrix(1,nj,1,nj);

	//set pollution levels, damages and circulation taxes per car to zero initially
	for (j = 1; j <= nj; j++) {
		co[j] = 0.0;
		hc[j] = 0.0;
		nox[j] = 0.0;
        totdamage[j] = 0.0;
        lctotdamage[j] = 0.0;
		damagepercar[j] = 0.0;
		tau[j] = 0.0;
	}
 
}


void dda_objective(int nummvec, double *mvec, double *diffs) {
    int j, g;
    double jsum, gsum;
    // the outside good is implicitly index j=0, so this represents the J+G-1 parameters/constraints
    
    for (g = 1; g <= ng; g++) {
        dda[0][g] = 1.0 * mvec[nj+g] * dda_base[0][g];
        for (j = 1; j <= nj; j++) {
            dda[j][g] = mvec[j] * mvec[nj+g] * dda_base[j][g];
        }
    }
    
    //compare dda vehicle counts to model
    for (j = 1; j <= nj; j++) {
        jsum=0.0;
        for (g = 1; g <= ng; g++) {
            jsum += dda[j][g];
        }
        diffs[j] = jsum - dda_jtarg[j];
    }
    
    //compare dda group pop counts to base
    for (g = 1; g <= ng; g++) {
        gsum = 0.0;
        for (j = 0; j <= nj; j++) gsum += dda[j][g];
        diffs[nj+g] = gsum - dda_gtarg[g];
    }
    
    
}


void writeout(int t) {
	//writes output for time t
	int j, g, m, dim1, dim2, dim3, dim4;
	double totdr1, totdr2, totgas1, totgas2, totco, tothc, totnox, evbound1, evbound2, temp1, temp2, temp3, temp4, temp5, temp6;
	double profit_bymake_1[NMAKE+1], profit_bymake_2[NMAKE+1], avgmarkup1[NMAKE+1], avgmarkup2[NMAKE+1], sales1[NMAKE+1], sales2[NMAKE+1];
	double mpgadj1[nj+1], mpgadj2[nj+1];
    double mvecguess[nj+ng+1], mvecdiffs[nj+ng+1];

	if (runmode==1 & splitused==0) {
		for (j=1;j<=njnew;j++) {
			mpgadj1[j]=mpg1[j]; 
			mpgadj2[j]=mpg2[j];
		}
		for (j=njnew+1;j<=nj;j++) {
			if (d1[j]<=usedsup1[j]) {
				mpgadj1[j]=mpg1[j];
				mpgadj2[j]=(usedsup2[j]/d2[j])*mpg2[j] + ((d2[j]-usedsup2[j])/d2[j])*mpg1[j];
			}
			else {
				mpgadj2[j]=mpg2[j];
				mpgadj1[j]=(usedsup1[j]/d1[j])*mpg1[j] + ((d1[j]-usedsup1[j])/d1[j])*mpg2[j];
			}
		}
	}
	else {	
		for (j=1;j<=nj;j++) {
			mpgadj1[j]=mpg1[j]; 
			mpgadj2[j]=mpg2[j];
		}
	}


	//******************************************************************************************************************
	// CAR LEVEL OUTPUTS
	
    //move total operating cost (rental plus fuel plus registration fee) into a j-indexed variable for reporting:
	for (dim1=1; dim1<=19; dim1++) {
		for (dim2=1; dim2<=2; dim2++) {
			for (dim3=1; dim3<=2; dim3++) {
				for (dim4=1; dim4<=7; dim4++) {
					tapj1[((dim1 - 1)*(NTYPE * NSIZE * NMAKE) + (dim2 - 1)*(NSIZE * NMAKE) + (dim3 - 1)*NMAKE + dim4)] = tap_1_1[dim1][dim2][dim3][dim4];
					tapj2[((dim1 - 1)*(NTYPE * NSIZE * NMAKE) + (dim2 - 1)*(NSIZE * NMAKE) + (dim3 - 1)*NMAKE + dim4)] = tap_1_2[dim1][dim2][dim3][dim4];
				}
			}
		}
	}

	// pollution levels in metric tons
	for (j = 1; j <= nj; j++) {
		co[j] = phi_co[j] * (d1[j] + d2[j]) * vmt[j] / 1000000;
		hc[j] = phi_hc[j] * (d1[j] + d2[j]) * vmt[j] / 1000000;
		nox[j] = phi_nox[j] * (d1[j] + d2[j]) * vmt[j] / 1000000;
        totdamage[j] = damage_co * co[j] + damage_hc * hc[j] + damage_nox * nox[j];
        lctotdamage[j] = totdamage[j] + proddamage[j]*(d1[j] + d2[j]);
	}
    
    // de-aggregate vehicles by income (holding shares proportional to baseline time 0)
    // target vectors:
    
    for (j = 1; j <= nj; j++) {
        dda_jtarg[j] = d1[j] + d2[j];
    }
    
    // group sizes
    for (g=1; g <= ng; g++) {
        dda_gtarg[g] = 0.0;
        for (j = 0; j <= nj; j++) dda_gtarg[g] += dda_base[j][g];
    }
    
    for (j=1; j<=nj+ng+1; j++) {
        mvecguess[j]=1.0;
    }
    broydn2(mvecguess,nj+ng,&retvaldda,dda_objective);
    if (retvaldda!=0) printf("\nERROR in income deaggregation");
    dda_objective(nj+ng,mvecguess,mvecdiffs);

    
    // overall file header (includes both car-level and aggregate data
    if (t==1) {
        fprintf(mainoutput, "run,baserun,t,j,age,type,size,make,mpg1,mpg2,cost1,cost2,rj1,rj2,d1,d2,pj1,pj2,pscrap1,pscrap2,hj1,hj2,dcde1,vmt,milesremaining1, fracremain0,dqde1,tapj1,tapj2, baserj1, basepj1, basepscrap1, commonmpg, phi_co, phi_hc, phi_nox, phi_co_when_new, phi_hc_when_new, phi_nox_when_new, tailpipecost, damagepercar, lcdamagepercar, tau, proptax, cartot_co, cartot_hc, cartot_nox, cartot_damage,cartot_lcdamage,d1_when_new,valj1,payment1,dda1,dda2,dda3,dda4,dda5,dda6,dda7,dda8,dda9,cardatafile, runmode, agg_d1, agg_d2, agg_gasoline1, agg_gasoline2, M_total_1, M_total_2, profit1, profit2, profitrental1, profitrental2, profit_adj1, profit_adj2, baseprofit1, baseprofit2, baseprofitrental1, baseprofitrental2, util1, util2, evbound1, evbound2, ev1, ev2, sumfocglobal, sumedglobal, retvalused, retvalnew, boundsviolations, agg_co, agg_hc, agg_nox, taxrev_tot_1, taxrev_tot_2, proptaxrev1, proptaxrev2 ,cafeCtarget, cafeTtarget, pavtarget,ddao1,ddao2,ddao3,ddao4,ddao5,ddao6,ddao7,ddao8,ddao9");
        
        //make-level data headers:
        for (m=1;m<=NMAKE;m++) fprintf(mainoutput,", cafeC_m%d",m);
        for (m=1;m<=NMAKE;m++) fprintf(mainoutput,", cafelimC_m%d",m);
        for (m=1;m<=NMAKE;m++) fprintf(mainoutput,", cafeCstat_m%d",m);
        for (m=1;m<=NMAKE;m++) fprintf(mainoutput,", lambdaC_m%d",m);

        for (m=1;m<=NMAKE;m++) fprintf(mainoutput,", cafeT_m%d",m);
        for (m=1;m<=NMAKE;m++) fprintf(mainoutput,", cafelimT_m%d",m);
        for (m=1;m<=NMAKE;m++) fprintf(mainoutput,", cafeTstat_m%d",m);
        for (m=1;m<=NMAKE;m++) fprintf(mainoutput,", lambdaT_m%d",m);
        
        for (m=1;m<=NMAKE;m++) fprintf(mainoutput,", cafeP_m%d",m);
        for (m=1;m<=NMAKE;m++) fprintf(mainoutput,", cafelimP_m%d",m);
        for (m=1;m<=NMAKE;m++) fprintf(mainoutput,", cafePstat_m%d",m);
        for (m=1;m<=NMAKE;m++) fprintf(mainoutput,", lambdaP_m%d",m);

        for (m=1;m<=NMAKE;m++) fprintf(mainoutput,", profit_bymake_1_m%d",m);
        for (m=1;m<=NMAKE;m++) fprintf(mainoutput,", sales1_m%d",m);
        for (m=1;m<=NMAKE;m++) fprintf(mainoutput,", avgmarkup1_m%d",m);
    }
    
    
    
    
	for (j = 1; j <= nj; j++) {
		fprintf(mainoutput, "\n\"%s\",\"%s\", %d, %d, %d, %d, %d, %d, %20.12lf, %20.12lf, %20.12lf, %20.12lf, %20.12lf, %20.12lf, %20.12lf, %20.12lf, %20.12lf, %20.12lf, %20.12lf, %20.12lf, %20.12lf, %20.12lf, %20.12lf, %20.12lf, %20.12lf, %20.12lf, %20.12lf, %20.12lf, %20.12lf, %20.12lf, %20.12lf, %20.12lf, %20.12lf, %20.12lf, %20.12lf, %20.12lf, %20.12lf, %20.12lf, %20.12lf, %20.12lf, %20.12lf, %20.12lf, %20.12lf, %20.12lf, %20.12lf, %20.12lf, %20.12lf, %20.12lf, %20.12lf, %20.12lf, %20.12lf, %20.12lf, %20.12lf, %20.12lf, %20.12lf, %20.12lf, %20.12lf, %20.12lf, %20.12lf, %20.12lf, %20.12lf",
			runname,
            baserunname,
            t,
            j,
            age[j],
			type[j],
			size[j],
			make[j],
			mpgadj1[j],
			mpgadj2[j],
			cost1[j],
			cost2[j],
			rj1[j],
			rj2[j],
			d1[j],
			d2[j],
			pj1[j],
			pj2[j],
			pscrap1[j],
			pscrap2[j],
            hj1[j],
            hj2[j],
			dcde1[j],
            vmt[j],
            milesremaining1[j],
            fracremain0[j],
            dqde1[j][j],
			tapj1[j],
			tapj2[j],
			baserj1[j],
			basepj1[j],
            basepscrap1[j],
			commonmpg[j],
            phi_co[j],
            phi_hc[j],
            phi_nox[j],
            phi_co_when_new[j],
            phi_hc_when_new[j],
            phi_nox_when_new[j],
            tailpipecost[j],
            damagepercar[j],
            lcdamagepercar[j],
            tau[j],
            proptax[j],
			co[j],
			hc[j],
			nox[j],
            totdamage[j],
            lctotdamage[j],
            d1_when_new[j],
            valj1[j],
            payment1[j][t],
            dda[j][1],
            dda[j][2],
            dda[j][3],
            dda[j][4],
            dda[j][5],
            dda[j][6],
            dda[j][7],
            dda[j][8],
            dda[j][9]
		);
	}
	
    
    
	//******************************************************************************************************************
	// AGGREGATE OUTPUTS
	
	//summary demand
	totdr1 = 0.0;
	totdr2 = 0.0;
	for (j=1;j<=nj;j++) {
		totdr1 += d1[j];
		totdr2 += d2[j];
	}

	//calculate total gasoline use
	totgas1 = 0;
	totgas2 = 0;
	for (j=1;j<=nj;j++) {
		totgas1 += d1[j] * vmt[j] / mpgadj1[j];
		totgas2 += d2[j] * vmt[j] / mpgadj2[j];
	}

	//calculate profits and average markup by firm
	for (j=1;j<=NMAKE;j++) {
		profit_bymake_1[j] = 0;
		profit_bymake_2[j] = 0;
		sales1[j] = 0;
		sales2[j] = 0;
	}
	for (j=1;j<=njnew;j++) {
		profit_bymake_1[make[j]] += (pj1[j] - cost1[j]) * d1[j];
		profit_bymake_2[make[j]] += (pj2[j] - cost2[j]) * d2[j];
		sales1[make[j]] += d1[j];
		sales2[make[j]] += d2[j];
    }
	for (j=1;j<=NMAKE;j++) {
		avgmarkup1[j] = profit_bymake_1[j] / sales1[j];
		avgmarkup2[j] = profit_bymake_2[j] / sales2[j];
	}

	//calculate total pollution
	totco = 0;
	tothc = 0;
	totnox = 0;
	for (j = 1; j <= nj; j++) {
		totco += co[j];
		tothc += hc[j];
		totnox += nox[j];
	}

	//calculate indirect utility
	indirectutil(dem_comp_5_1[1],dem_comp_5_1[2],dem_comp_5_2[1],dem_comp_5_2[2]);

    //calculate bound for EV
    temp1 = 0;
    temp2 = 0;
    for (j=1;j<=nj;j++) {
        temp1 -= (based1[j] * (rj1[j] - baserj1[j]) + 0.5 * (rj1[j] - baserj1[j]) * (d1[j] - based1[j]));
        temp2 -= (based2[j] * (rj2[j] - baserj2[j]) + 0.5 * (rj2[j] - baserj2[j]) * (d2[j] - based2[j]));
    }
    evbound1 = temp1 + profit_adj1;
    evbound2 = temp2 + profit_adj2;
    

    
    //note: aggregates go to the right of the car data -- the long series of commas skips all those columns
    fprintf(mainoutput, "\n\"%s\",\"%s\",%d,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,\"%s\",%d,%20.12lf,%20.12lf,%20.12lf,%20.12lf,%20.12lf,%20.12lf,%20.12lf,%20.12lf,%20.12lf,%20.12lf,%20.12lf,%20.12lf,%20.12lf,%20.12lf,%20.12lf,%20.12lf,%20.12lf,%20.12lf,%20.12lf,%20.12lf,%20.12lf,%20.12lf,%20.12lf,%20.12lf,%d,%d,%d,%20.12lf,%20.12lf,%20.12lf,%20.12lf,%20.12lf,%20.12lf,%20.12lf,%20.12lf,%20.12lf,%20.12lf,%20.12lf,%20.12lf,%20.12lf,%20.12lf,%20.12lf,%20.12lf,%20.12lf,%20.12lf,%20.12lf",
            runname,
            baserunname,
            t,
            
            cardatafile,
            runmode,
            totdr1,
            totdr2,
            totgas1,
            totgas2,
            M_total_1,
            M_total_2,
            profit1,
            profit2,
            profitrental1,
            profitrental2,
            profit_adj1,
            profit_adj2,
            baseprofit1,
            baseprofit2,
            baseprofitrental1,
            baseprofitrental2,
            util1[t][runmode],
            util2[t][runmode],
            evbound1,
            evbound2,
            ev1[t],
            ev2[t],
            sumfocglobal,
            sumedglobal,
            retvalused,
            retvalnew,
            boundsviolations,
            totco,
            tothc,
            totnox,
            taxrev_tot_1,
            taxrev_tot_2,
            proptaxrev1,
            proptaxrev2,
            cafeCtarget[t],
            cafeTtarget[t],
            pavtarget[t],
            dda[0][1],
            dda[0][2],
            dda[0][3],
            dda[0][4],
            dda[0][5],
            dda[0][6],
            dda[0][7],
            dda[0][8],
            dda[0][9]
    );
    
    //append make-level data to the output row

    for (m=1;m<=NMAKE;m++) fprintf(mainoutput,",%20.12lf",cafeC[m]);
    for (m=1;m<=NMAKE;m++) fprintf(mainoutput,",%20.12lf",cafelimC[m]);
    for (m=1;m<=NMAKE;m++) fprintf(mainoutput,",%20.12lf",cafeCstat[m]);
    for (m=1;m<=NMAKE;m++) fprintf(mainoutput,",%20.12lf",lambdaC[m]);
    
    for (m=1;m<=NMAKE;m++) fprintf(mainoutput,",%20.12lf",cafeT[m]);
    for (m=1;m<=NMAKE;m++) fprintf(mainoutput,",%20.12lf",cafelimT[m]);
    for (m=1;m<=NMAKE;m++) fprintf(mainoutput,",%20.12lf",cafeTstat[m]);
    for (m=1;m<=NMAKE;m++) fprintf(mainoutput,",%20.12lf",lambdaT[m]);

    for (m=1;m<=NMAKE;m++) fprintf(mainoutput,",%20.12lf",cafeP[m]);
    for (m=1;m<=NMAKE;m++) fprintf(mainoutput,",%20.12lf",cafelimP[m]);
    for (m=1;m<=NMAKE;m++) fprintf(mainoutput,",%20.12lf",cafePstat[m]);
    for (m=1;m<=NMAKE;m++) fprintf(mainoutput,",%20.12lf",lambdaP[m]);

    for (m=1;m<=NMAKE;m++) fprintf(mainoutput,",%20.12lf",profit_bymake_1[m]);
    for (m=1;m<=NMAKE;m++) fprintf(mainoutput,",%20.12lf",sales1[m]);
    for (m=1;m<=NMAKE;m++) fprintf(mainoutput,",%20.12lf",avgmarkup1[m]);
}


void advance(int t) {
	int j;
    double temppaytotal1[nj+1],temppaytotal2[nj+1];
	
    
    // car payment model
    for (j=1; j<=njnew; j++) {
        d1_when_new[j] = d1[j];  //store demand when new for the current new cars
        d2_when_new[j] = d2[j];
    }
  
    // reconcile any differences in payments and purchase price (e.g. coming from endogeneity of scrap rates and mpg adjustments)
    for (j=1; j<=nj; j++) {
        if (age[j]==0) {
            temppaytotal1[j] = payment1[j][t+age[j]];
            temppaytotal2[j] = payment2[j][t+age[j]];
        }
        else {
            temppaytotal1[j-njnew*age[j]] += payment1[j][t+age[j]] / pow(1.0+drate,age[j]);
            temppaytotal2[j-njnew*age[j]] += payment2[j][t+age[j]] / pow(1.0+drate,age[j]);
        }
    }
    
    for (j=1; j<=nj; j++) {
        payment1[j][t+age[j]] = payment1[j][t+age[j]] * pj1[j-njnew*age[j]]/temppaytotal1[j-njnew*age[j]] ;
        payment2[j][t+age[j]] = payment2[j][t+age[j]] * pj2[j-njnew*age[j]]/temppaytotal2[j-njnew*age[j]] ;
    }
    
    
	// used car attributes to move with the car through time: maxsup, mpg, income
	for (j=nj;j>=njnew+1;j--) {
		mpg1[j] = mpg1[j-njnew];
		mpg2[j] = mpg2[j-njnew];
		maxsup1[j] = d1[j-njnew];
		maxsup2[j] = d2[j-njnew];
        d1_when_new[j] = d1_when_new[j-njnew];
        d2_when_new[j] = d2_when_new[j-njnew];
    }
    
    // store values from time t-1
    for (j=1; j<=nj ;j++) {
        lastpj1[j]=pj1[j];
        lastpj2[j]=pj2[j];
        lastd1[j] = d1[j];
        lastd2[j] = d2[j];
        lastexppj1[j] = exppj1[j];
        lastexppj2[j] = exppj2[j];
        lastexppscrap1[j] = exppscrap1[j];
        lastexppscrap2[j] = exppscrap2[j];
        lastexphj1[j] = exphj1[j];
        lastexphj2[j] = exphj2[j];
    }
    

	M_1 = M_1 * (1.0+M_growth);
	M_2 = M_2 * (1.0+M_growth);

	//exog tech change in basempg of cars (free gain in mpg each time step, cost function shifts right)
    //and in production damage
	for (j=1;j<=njnew;j++) {
		basempg1[j] += exogmpgtech * basempg1[j];
		basempg2[j] += exogmpgtech * basempg2[j];
        proddamage[j] = proddamage[j]*0.9193;
	}
    
    
    //pull emissions factors and tailpipe cost in from the j x t matrix (calculated in definecar)
    for (j = 1; j <= nj; j++) {
        phi_co[j] = jt_phi_co[j][t+1];
        phi_hc[j] = jt_phi_hc[j][t+1];
        phi_nox[j] = jt_phi_nox[j][t+1];
        phi_co_when_new[j] = jt_phi_co_when_new[j][t+1];
        phi_hc_when_new[j] = jt_phi_hc_when_new[j][t+1];
        phi_nox_when_new[j] = jt_phi_nox_when_new[j][t+1];
    }
    for (j = 1; j <= njnew; j++) {
        tailpipecost[j] = jt_tailpipecost[j][t+1];
    }
    
    // variables to reinitialize to zero: demands
    for (j=1;j<=nj;j++) {
        d1[j] = 0;
        d2[j] = 0;
    }
}
